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1. INTRODUCTION 


As part of the continuing effort at NASA/Lewis to improve both the durability and 
reliability of hot section Earth- to- Orbit engine components, significant enhancements must 
be made in existing finite element and finite difference methods, and advanced techniques, 
such as the boundary element method, must be explored. Despite this considerable effort, 
the accurate determination of transient thermal stresses in these hot section components 
remains one of the most difficult problems facing engine design/analysts. For these prob- 
lems, the temperature distribution is strongly influenced by the external hot gas flow, 
the internal cooling system, and the structural deformation. Currently, experimentally- 
determined film coefficients and ambient temperatures are required for use as boundary 
conditions for the thermal stress analysis of the structural component. The determina- 
tion of these coefficients is obviously an expensive and time-consuming task. Recently an 
attempt was made by Gladden (1989) to use a finite difference-based Navier-Stokes code 
to approximate the thermal boundary conditions, and to then input these into a finite 
element structural analysis package. However, the most effective way to deal with this 
problem is to develop a completely integrated solid mechanics, fluid mechanics, and heat 
transfer approach. 

In the present work, the boundary element method (BEM) is chosen as the basic 
analysis tool principally because the critical surface variables (i.e., temperature, flux, dis- 
placement, traction) can be very precisely determined with a boundary-based discretization 
scheme. Additionally, model preparation is considerably simplified compared to the more 
fam iliar domain-based methods. Furthermore, the hyperbolic character of high speed flow 
is captured through the use of an analytical fundamental solution, eliminating the depen- 
dence of the solution on the discretization pattern. The price that must be paid in order 
to realize these advantages is that any BEM formulation requires a considerable amount 
of analytical work, which is typically absent in the other numerical methods. 


1 



This report details all of the research accomplishments of a multi-year program, com- 
mencing in March 1986, aimed toward the development of a boundary element formulation 
for the study of hot fluid-structure interaction in Earth-to-Orbit engine hot section com- 
ponents, It should be noted that this work represents approximately four man-years of 
funding from NASA/Lewis. Most of that effort expended under this program has been di- 
rected toward the examination of fluid flow, since boundary element methods for fluids are 
at a much less developed state. Recently, however, significant strides have been made, not 
only in the analysis of thermoviscous fluids, but also in the solution of the fluid-structure 
interaction problem. 

Early in the research program, a two-dimensional boundary element formulation was 
developed for the time- dependent response of a thermoelastic solid. This effort resulted 
in the first time domain, boundary-only implementation for this class of problems. Since 
volume discretization is completely eliminated and surface transient thermal stresses can 
be captured very accurately, the new approach provides distinct advantages over standard 
finite element methods. 

Meanwhile, the initial fluid formulations that were developed, based upon Stokes fun- 
damental solutions, provided solutions in the low-to-moderate Reynolds number range. 
For creeping flow, these reduce to boundary-only techniques. As the fluid velocities are in- 
creased, volume discretization is required, however the solutions are typically very precise, 
particularly in the determination of surface quantities. At very high speed, these formu- 
lations are less effective, because the Stokes fundamental solutions no longer embody the 
character of the flow field which becomes dominated by convection. 

This led to the development of convective viscous integral formulations based upon Os- 
een fundamental solutions. Since the new convective kernel functions, that were developed 
as a part of this effort, contain more of the physics of the problem, boundary element so- 
lutions can now be obtained at very high Reynolds number. Flow around obstacles can be 
solved approximately with an efficient linearized boundary-only analysis or more exactly 
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by including all of the nonlinearities present in the neighborhood of the obstacle. This 
perhaps represents the major accomplishment of the present program. 

The other significant development has been the creation of a comprehensive fluid- 
structure interaction capability within a boundary element computer code. This new 
facility is implemented in a completely general manner, so that quite arbitrary geometry, 
material properties and boundary conditions may be specified. Thus, a single analysis 
code can be used to run structures-only problems, fluids-only problems, or the combined 
fluid-structure problem. In all three cases, steady or transient conditions can be selected, 
with or without thermal effects. Nonlinear analyses can be solved via direct iteration or by 
employing a modified Newton-Raphson approach. However, it should be emphasized that 
the existing program is primarily a research code. Significant additional effort is needed 
to develop a practical engineering analysis tool. 

In the next section, a brief review of the recent applicable boundary element literature 
is presented. This is followed by the development of integral formulations for the ther- 
moelastic solid in Section 3 and for the thermoviscous fluid in Section 4. A number of 
detailed numerical examples are included at the end of these two sections to validate the 
formulations and to emphasize both the accuracy and generality of the implementation. 
Then, in Section 5, the fluid-structure interaction facility is discussed. Once again, several 
examples axe provided to highlight this unique capability. It should be noted that all of 
the results presented in this report were run on a desktop SUN SPARCstation 1. Section 6 
contains a collection of potential boundary element applications that have been uncovered 
as a result of work related to the present grant. For most of those problems, satisfactory 
analysis techniques do not currently exist. The remaining sections summarize the progress 
achieved to date, and specify the future direction. Tables and figures appear at the end of 
each section, while references are provided in Appendix A. 
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2. LITERATURE REVIEW 

Very little has appeared in the literature on the analysis of coupled thermoviscous fluid- 
structure problems via the boundary element method. However, a number of publications 
have addressed the fluid and structure separately. 

In general, the solid portion of the problem has been addressed to a much greater 
degree. For example, a boundary-only steady-state thermoelastic formulation was initially 
presented by Cruse et al (1977) and Rizzo and Shippy (1977). Recently, the present 
authors developed and implemented the quasistatic counterpart (Dargush, 1987, Dargush 
and Banerjee, 1989b, 1990a, 1990b), which is presented in detail in Section 3. Others, 
notably Sharp and Crouch (1986) and Chaudouet (1987), introduce volume integrals, to 
represent the equivalent thermal body forces. A similar domain based approach was taken 
earlier by Banerjee and Butterfield (1981) in the context of the analogous geomechanical 
problem. 

An extensive review of the applications of integral formulations to viscous flow prob- 
lems was included in a previous annual report (Dargush et al, 1987), and will not be 
repeated here. Interestingly, only a few groups of researchers are actively pursuing the 
further development of boundary elements for the analysis of viscous fluids. The work re- 
ported in Piva and Morino (1987) and Piva et al (1987) focuses heavily on the development 
of fundamental solutions and integral formulations with little emphasis on implementation. 
On the other hand, Tosaka and Kakuda (1986, 1987), Tosaka and Onishi (1986) have im- 
plemented single region boundary element formulations using approximate incompressible 
fundamental solutions. This latter group has developed sophisticated non-linear solution 
algorithms, and consequently, are able to demonstrate moderately high Reynolds num- 
ber solutions. Meanwhile, Dargush and Banerjee (1991a, 1991b) present general purpose 
steady and time- dependent boundary element methods for moderate Reynolds number 
flows. 
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The most recent work from the above researchers has been collected into a volume en- 
titled Developments in BEM - Volume 6: Nonlinear Probl ems of Fluid Dynamics, edited 
by Banerjee and Morino. Contributions from Wu and Wang, and Bush and Tanner are also 
included, along with two chapters from the present co-authors. The volume, published by 
Elsevier Applied Science Publishers became available in mid-1990, and provides a state- 
of-the-art review of boundary element fluid dynamics. However, it should be noted that 
the convective thermoviscous formulations of Section 4 are not included. These represent 
a significant further advancement which permit solutions for high Reynolds number flows. 
Interestingly, the basis for much of this latter development is actually work done early in 
this century by Oseen (1911, 1927). 

For analysis of the interaction problem, a boundary element thermoelastic solid repre- 
sentation must be coupled with a suitable thermoviscous fluid formulation. Only Dargush 
and Banerjee (1988,1989a) have tackled this problem. These two papers provide a sum- 
mary of the early work performed under this grant. 
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3. INTEGRAL FORMULATION FOR SOLIDS 

3.1 Introduction 

In the current section, a surface only time domain boundary element method (BEM) 
will be described for a thermoelastic body under quasistatic loading. Thus, transient heat 
conduction is included, but inertial effects are ignored. This BEM was first developed as 
part of the work performed during the second year (1987) of this grant. Since that time a 
number of improvements and extensions have been incorporated. During 1989, the algo- 
rithms for numerical integration have been made more efficient as well as more accurate, 
and a comprehensive PATRAN interface has been added to aid in the post-processing of 
the boundary element results. Additionally, a streamlined approach for uncoupled ther- 
moelasticity was introduced (Dargush and Banerjee, 1989b). In 1990, boundary elements 
with a quartic variation of the field variables were implemented. These elements are par- 
ticularly well suited for problems involving the bending of components (Deb and Banerjee, 
1989). 

Details of the integral formulation for 2D plane strain is presented below. (Problems 
of plane stress can be handled via a simple change in material parameters.) Separate sub- 
sections present the governing differential equations, the integral equations, an overview 
of the numerical implementation, and a couple of simple examples. Similar formulations 
have also been developed for three-dimensional (Dargush and Banerjee, 1990a) and ax- 
isymmetric problems (Dargush and Banerjee, 1990b). 


3.2 Governing Equations 

With the solid assumed to be a linear thermoelastic medium, the governing differential 


equations for transient thermoelasticity can be written 


(A + - (W + 2 «>«lr = o 


dxidxj 


dxjdxj 

d$_ 
dt 


dx{ 


eo , e 2 e 

pc c — = k 


dxjdxj 


( 3 . 1 a) 


( 3 . 16 ) 
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where 


ui displacement vector 
9 temperature 
t time 

x% Lagrangian coordinate 
k thermal conductivity 
p mass density 

c e specific heat at constant deformation 
A, p Lame constants 
a coefficient of thermal expansion 

Standard indicia! notation has been employed with summations indicated by repeated 
indices. For two-dimensional problems considered herein, the Latin indices i and j vary 
from one to two. 

Note that (3.1b) is the energy equation and that (3.1a) represents the momentum 
balance in terms of displacements and temperature. The theory portrayed by the above 
set of equations, formally labeled uncoupled quasistatic thermoelasticity, can be derived 
from thermodynamic principles. (See Boley and Weiner (1960) for details.) In developing 
(3.1), the dynamics effects of interia have been ignored. 


3.3 Integral Representations 

Utilizing equation (3.1) for the solid along with a generalized form of the reciprocal 
theorem, permits one to develop the following boundary integral equation: 


C(3a{£) u p{£ j 0 


l 


9f)a*t0(X,t) - fp a * Up(X,t) 


dS(X). 


(3.2) 


where 

a,/3 indices varying from 1 to 3 
s surface of solid 


7 


u Q ,t a generalized displacement and traction 

tic, = [«1 «2 ^] T 

to = [<1 <2 ?] T 
6,q temperature, heat flux 

gap, fop generalized displacement and traction kernels 

c aP , constants determined by the relative smoothness of s at £ 

and, for example 

g a p*ta = / 9ap{ x >t'i£i T )ta( x > T )d T 

JO 

denotes a Riemann convolution integral. The kernel functions g Q p and f Q p are derived from 
the fundamental infinite space solutions of (3.1). 

In principle, at each instant of time progressing from time zero, this equation can be 
written at every point on the boundary. The collection of the resulting equations could then 
be solved simultaneously, producing exact values for all the unknown boundary quantities. 
In reality, of course, discretization is needed to limit this process to a finite number of 
equations and unknowns. Techniques useful for the discretization of (3.2) are the subject 
of the following section. 

3.4 Numerical Implementation 

3.4.1 Introduction 

The boundary integral equation (3.2), developed in the last section, is an exact state- 
ment. No approximations have been introduced other than those used to formulate the 
boundary value problem. However, in order to apply (3.2) for the solution of practical en- 
gineering problems, approximations are required in both time and space. In this section, 
an overview of a general-purpose, state-of-the-art numerical implementation is presented. 
Many of the features and techniques to be discussed, in this section, were developed previ- 
ously for elastostatics (e.g., Banerjee et al, 1985, 1988), and elastodynamics (e.g., Banerjee 
et al, 1986; Ahmad and Banerjee, 1988), but are here adapted for thermoelastic analysis. 


8 



3.4.2 Temporal Discretization 

Consider, first, the time integrals represented in (3.2) as convolutions. Clearly, without 
any loss of precision, the time interval from zero to t can be divided into N equal increments 
of duration At. 

By assuming that the primary field variables, tp and up, are constant within each At 
time increment, these quantities can be brought outside of the time integral. That is, 

N ynAt 

gpa*tp(X,i) = J2t%{X) gpa{X -Z,t-T)dT (3.3a) 

J(n- 1)A* 


N AnAt 

ffic * Ufi(X, t)=T U%(X) / fpa(X -t,t- r)dr (3.36) 

where the superscript on the generalized tractions and displacements, obviously, represents 
the time increment number. Notice, also, that, within an increment, these primary field 
variables are now functions of position only. Next, since the integrands remaining in 
(3.3) are known in explicit form from the fundamental solutions, the required temporal 
integration can be performed analytically, and written as 

pnAt 

G»+ 1 ~ n (X-0= 9p a (X-t,t-T)dT (3.4a) 

J (n— l)At 


pnAt 

F" +1 -"(X-0=/ ffia(X-U-r)dr. (3.46) 

J (n— l)At 


These kernel functions, G% a (X~Z) and F$ Q {X -{), are detailed in Appendix B.l. Combining 
(3.3) and (3.4) with (3.2) produces 


c*,(0«?(0 = £ / te +1 -"(^ - otp(x) - F” 0 +1 ~ n (x - Oujftx) 

n=l Js L 


dS(X) f 


(3.5) 


which is the boundary integral statement after the application of the temporal discretiza- 
tion. 
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3.4.3 Spatial Discretization 

With the use of generalized primary variables and the incorporation of a piecewise 
constant time stepping algorithm, the boundary integral equation (3.5) begins to show 
a strong resemblance to that of elastostatics, particularly for the initial time step (i.e. , 
N = l). In this subsection, those similarities will be exploited to develop the spatial 
discretization for the uncoupled quasistatic problem with two-dimensional geometry. This 
approximate spatial representation will, subsequently, permit numerical evaluation of the 
surface integrals appearing in (3.5). The techniques described here, actually, originated in 
the finite element literature, but were later applied to boundary elements by Lachat and 
Watson (1976). 

The process begins by subdividing the entire surface of the body into individual ele- 
ments of relatively simple shape. The geometry of each element is, then, completely defined 
by the coordinates of the nodal points and associated interpolation functions. That is, 

X(() = x i(C) = N w (()Ziw (3-6) 

with 

( intrinsic coordinates 
N w shape functions 
n w nodal coordinates 

and where w is an integer varying from one to W , the number of geometric nodes in the 
element. Next, the same type of representation is used, within the element, to describe 
the primary variables. Thus, 

«S(C) = AUCKu, ( 3 ' 7fl ) 

a o = auoCw ( 3 - 7fc ) 

in which and are the nodal values of the generalized displacement and tractions, 
respectively, for time step n. Also, in (3.7), the integer w varies from one to Q, the total 
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number of functional nodes in the element. From the above, note that the same number 
of nodes, and consequently shape functions, are not necessarily used to describe both the 
geometric and functional variations. Specifically, in the present work, the geometry is 
exclusively defined by quadratic shape functions. In two-dimensions, this requires the use 
of three-noded line elements. On the other hand, the variation of the primary quantities can 
be described, within an element, by linear, quadratic or quartic shape functions. For each 
quartic element, two additional quarter-point nodes are automatically generated by the 
program. It should be noted that the introduction of quartic elements this past year, also 
provides the foundation for the development of a p-adaptive boundary element capability. 

Once the spatial discretization has been accomplished and the body has been subdi- 
vided into M elements, the boundary integral equation can be rewritten as 

n , m , r 

e*,(O«F(0 = £ £ / - ONUC)t n p. 

n=l ^ m= 1 JS ™ L 

- F” +1 ~ n (X ( 0 - 0*U0«&, dS(X(C))| (3-8) 

where 

M 

S = 2 Sm- 

m= 1 

In the above equation, t n 0w and are nodal quantities which can be brought outside the 
surface integrals. Thus, 

C „»( »<«) = E { E / CJJ'-WO - OWWB) 

n= 1 ^ m=l JSm 

- J s F^ +1 - n (X(0 - QNUQdSWQ)} (3.9) 

The positioning of the nodal primary variables outside the integrals is, of course, a key 
step since now the integrands contain only known functions. However, before discussing 
the techniques used to numerically evaluate these integrals, a brief discussion of the sin- 
gularities present in the kernels and F£ a is in order. 

The fundamental solutions to the uncoupled quasistatic problem contain singularities 
when the load point and field point coincide, that is, is when r = 0. The same is true of G 0a 
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and Fp Q , since these kernels axe derived directly from the fundamental solutions. Series 
expansions of terms present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A number of observations concerning the results of these expansions should be men- 
tioned. First, as would be expected F*p has a stronger level of singularity than does the 
corresponding G l Q(3} since an additional derivative is involved in obtaining F^p from G l a p . 
Second, the coupling terms do not have as a high degree of singularity as do the corre- 
sponding non-coupling terms. Third, all of the kernel functions for the first time step could 
actually be rewritten as a sum of steady-state and transient components. That is, 

G 1 _ss /-* i tr s~il 
a(3 — + U a{3 

Fi P =“ s F a p + tr F l ap . 

Then, the singularity is completely contained in the steady-state portion. Furthermore, 
the singularity in Gjj and F- is precisely equal to that for elastostatics, while Ggg and Fgg 
singularities are identical to those for potential flow. (For two-dimensions, the subscript 
6 equals three.) This observation is critical in the numerical integration of the F a p kernel 
to be discussed in the next subsection. However, from a physical standpoint, this means 
that, at any time t, the nearer one moves toward the load point, the closer the quasistatic 
response field corresponds with a steady-state field. Eventually, when the sampling and 
load points coincide, the quasistatic and steady-state responses are indistinguishable. As 
a final item, after careful examination of Appendix B.l, it is evident that the steady-state 
components in the kernels G n op and F^,, with n > 1, vanish. In that case, all that remains 
is a transient portion that contains no singularities. Thus, all singularities reside in the 
“G aP and “F a0 components of G l a0 and F^, respectively. 

3.4.4 Numerical Integration 

Having clarified the potential singularities present in the coupled kernels, it is now 
possible to consider the evaluation of the integrals in equation (3.9). That is, for any 
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element m, the integrals 



g at+1 -"(*«) -0MC)dS{x(0) 



F N+1 - n (X(0-t)NUOdS(X(0) 


(3.10a) 

(3.106) 


will be examined. To assist in this endeavor, the following three distinct categories can be 


identified. 


(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but only non-singular or weakly singular integrals 

axe involved. 

(3) The point £ lies on the element m, and the integral is strongly singular. 

In practical problems involving many elements, it is evident that most of the integration 
occurring in equation (3.9) will be of the category (1) variety. In this case, the integrand is 
always non-singular, and standard Gaussian quadrature formulas can be employed. Sophis- 
ticated error control routines are needed, however, to minimize the computational effort 
for a certain level of accuracy. This non-singular integration is the most expensive part of 
a boundary element analysis, and, consequently, must be optimized to achieve an efficient 
solution. In the present implementation, error estimates, based upon the work of Stroud 
and Secrest (1966), are employed to automatically select the proper order of the quadrature 
rule. Additionally, to improve accuracy in a cost-effective manner, a graded subdivision 
of the element is incorporated, especially when £ is nearby. For two-dimensional prob- 
lems, the integration order varies from two to twelve, within each of up to four element 
subdivisions. 

Turning next to category (2), one finds that again Gaussian quadrature is applicable, 
however, a somewhat modified scheme must be utilized to evaluate the weakly singular 
integrals. This is accomplished in two-dimensional elements via suitable subsegmentation 
along the length of the element so that the product of shape function, Jacobian and kernel 
remains well behaved. 
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Unfortunately, the remaining strongly singular integrals of category (3) exist only in 
the Cauchy principal value sense and cannot, in general, be evaluated numerically, with 
sufficient precision. It should be noted that this apparent stumbling block is limited to the 
strongly singular portions, “Fy and ••Fee, of the F l ap kernel. The remainder of F^ p , including 
tr Ffj and tr Fg e , can be computed using the procedures outlined for category (2). However, 
as will be discussed in the next subsection, even category (3) **Fy and “Fee kernels can be 
accurately determined by employing an indirect ‘rigid body’ method originally developed 
by Cruse (1974). 


3.4.5 Assembly 

The complete discretization of the boundary integral equation, in both time and space, 
has been described, along with the techniques required for numerical integration of the 
kernels. Now, a system of algebraic equations can be developed to permit the approximate 
solution of the original quasistatic problem. This is accomplished by systematically writing 
(3.9) at each global boundary node. The ensuing nodal collocation process, then, produces 
a global set of equations of the form 


N 

£ 

n — 1 


qN - f 1 — n 




jpiiV-fl — 


"]{«"}) = {o}, 


(3.11) 


where 

[£jv+i-n] unassembled matrix of size ( d + 1 )P x (d + 1)Q, with coefficients determined 
from (3.10a) 

[F N+ i-nj assembled matrix of size (d+\)P x (d+l)P, with coefficients determined from 
(3.10b) and c Pa included in the diagonal blocks 

{t n } global generalized nodal traction vector with ( d + 1)Q components 

{u n } global generalized nodal displacement vector with ( d+ 1 )P components 


{0} null vector with (d+ 1 )P components 
P total number of global functional nodes 
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Q — 5Zm=l Am 


A m number of functional nodes in element m 
d dimensionality of the problem. 

In the above, recall that the terms generalized displacement and traction refer to the 
inclusion of the temperature and flux, respectively, as the (cf + 1) component at any point. 
Consider, now, the first step. Thus, for N = 1, equation (3.11) becomes 

[G 1 ]^ 1 } — [i rl ]{u 1 > = {0}. (3.12) 

However, at this point the diagonal block of [F 1 ] has not been completely determined due to 
the strongly singular nature of 3S Fij and 33 Fee- Following Cruse (1974) and, later, Banerjee 
et al (1986) in elastodynamics, these diagonal contributions can be calculated indirectly 
by imposing a uniform ‘rigid body’ generalized displacement field on the same body, but 
under steady-state conditions. Then, obviously, the generalized tractions must be zero, 
and 

["F]{1} = {0}, ( 3 - 13 ) 

where {1} is a vector symbolizing a unit uniform motion. Using (3.13), the desired diagonal 
blocks, " Fij and 33 Fee, can be obtained from the summation of the off-diagonal terms of 
[ J *F]. The remaining transient portion of the diagonal block is non-singular, and hence 
can be evaluated to any desired precision. After summing the steady-state and transient 
contributions, (3.12) is once again written as 

[G 1 ]{t 1 }-[F 1 ]{a 1 } = {0} ) (3.14) 

but now the evaluation of [F 1 ] is complete. 

In a well-posed problem, at time At, the set of global generalized nodal displacements 
and tractions will contain exactly ( d + 1 )P unknown components. Then, as the final stage 
in the assembly process, equation (3.14) can be rearranged to form 

[A l ]{x l ] = miy 1 }, (3.15) 


15 


in which 


{x 1 } unknown components of {u 1 } and {i 1 } 

{y 1 } known components of {u 1 } and {t 1 } 

[A 1 ], [B 1 ] associated matrices 

3.4.6 Solution 

To obtain a solution of (3.15) for the unknown nodal quantities, a decomposition 
of matrix [A 1 ] is required. In general, [A 1 ] is a densely populated, unsymmetric matrix. 
The out-of-core solver, utilized here, was developed originally for elastostatics from the 
LINPACK software package (Dongarra et al, 1979) and operates on a submatrix level. 
Within each submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [.A 1 ] is stored in a direct-access file 
for reuse in subsequent time steps. Backsubstitution then completes the determination of 
{x 1 }. Additional information on this solver is available in Banerjee et al (1985). 

After turning from the solver routines, the entire nodal response vectors, {u 1 } and 
{t 1 }, at time At are known. For solutions at later times, a simple marching algorithm is 
employed. Thus, from (3.11) with N = 2, 

[G l ]{P} - [F 1 ]^ 1 } + [G 1 ]^ 2 } - [F 1 ]^ 2 } = {0}. (3.16) 

Assuming that the same set of nodal components are unknown as in (3.14) for the first 
time step, equation (3.16) is reformulated as 

[yF]{x 2 } = [Id]{y 2 } - [G 2 ]{P} + [F 2 ]^ 1 }. (3.17) 

Since, at this point, the right-hand side contains only known quantities, (3.17) can be 
solved for {x 2 }. However, the decomposed form of [A 1 ] already exists on a direct-access file, 
so only the relatively inexpensive backsubstitution phase is required for the solution. 

The generalization of (3.17) to any time step N is simply 

[A 1 }{x N } = m{y"}-Y t 


f[G N+1_n ]{t n } - [F N+1-n ]{ti n }J 


(3.18) 


16 


in which the summation represents the effect of past events. By systematically storing 
all of the matrices and nodal response vectors computed during the marching process, 
surprisingly little computing time is required at each new time step. In fact, for any 
time step beyond the first, the only major computational task is the integration needed 
to form [G N ] and [F N ]. Even this process is somewhat simplified, since now the kernels 
are non-singular. As a result, reduced subsegmentation and gaussian integration order is 
appropriate. Also, as time marches on, the effect of events that occurred during the first 
time step diminishes. Consequently, the terms containing [G w ] and [ F N ] will eventually 
become insignificant compared to those associated with recent events. Once that point is 
reached, further integration is unnecessary, and a significant reduction in the computing 
effort per time step can be achieved. 

It should be emphasized that the entire boundary element method developed, in this 
section, has involved surface quantities exclusively. A complete solution to the well-posed 
linear uncoupled quasistatic problem, with homogeneous properties, can be obtained in 
terms of the nodal response vectors, without the need for any volume discretization. In 
many practical situations, however, additional information, such as, the temperature at 
interior locations or the stress at points on the boundary, is required. The next subsection 
discusses the calculations of these quantities. 

3.4.7 Interior Quantities 

Once equation (3.18) is solved, at any time step, the complete set of primary nodal 
quantities, {u n } and {< N }, is known. Subsequently, the response at points within the body 
ran be calculated in a straightforward manner. For any point £ in the interior, the gener- 
alized displacement can be determined from (3.9) with cp a = 6p a . That is, 

jv r m r r 

«s(o = E £ kJ g" + 1 ~"(*( o - ojmo^wo) 

n=l ^ m - 1 *• • >Sm 

- «2L, j s *& +, “"WC) - OAUC)dS(X(C))] }• (3.19) 

Now, all the nodal variables on the right-hand side are known, and, as long as, t is not on 
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the boundary, the kernel functions in (3.19) remain non-singular. However, when £ is on the 
boundary, the strong singularity in ‘‘ a Fp a prohibits accurate evaluation of the generalized 
displacement via (3.19), and an alternate approach is required. The apparent dilemma is 
easily resolved by recalling that the variation of surface quantities is completely defined 
by the elemental shape functions. Thus, for boundary points, the desired relationship is 
simply 

*£«) = JU0«2L (3 ' 20) 


where N u { C) are the shape functions for the appropriate element and C are the intrinsic 
coordinates corresponding to £ within that element. Obviously, from (3.20), neither in- 
tegration nor the explicit contribution of past events are needed to evaluate generalized 
boundary displacements. 

In many problems, additional quantities, such a heat flux and stress, are also important. 


The boundary integral equation for heat flux, can be written 

n { m r r 

rf'io-E E w iff-wo-wwio) 

n=l ^ m=l *' 

- J s D% e f- n (X(0 - 0^(C)d5(X(0)] |. 


(3.21) 


where 


X(C)-0 = ~k 

W0-0 = -* 


dG" pe (X( 0-0 

% 

dFfajxj Q - Q 

% 


(3.21a) 

(3.216) 


This is valid for interior points, whereas, when £ is on the boundary, the shape functions 
can again be used. In this latter case, 


auo?" = * i(o<?r(o ( 3 - 22a ) 

= ( 3 - 22fc ) 

which can be solved for boundary flux. Meanwhile, interior stresses can be evaluated from 

N , M f r 

<(0 = £ £ KL / <f~"(*(C) - Z)NM)dS(X(0) 

n=l m=l L JS ™ 

~ - 0A r U0^(^(0)l } (3.23) 

J Sm J ' 
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in which 


EfiiiiM 0 - 0 


1 - 2v 


Utj ■ 


8U 



'8G% 

dti ) 

| - ftuCfa 

(3.23a) 

+d 

\ 

♦2) 

— / 36ij Fp & , 

(3.236) 


vpiiiHQ - 0 = Y^2u Sii ~dit 

with v representing the Poisson ratio and (5 = (3 A -f 2y)a. Equation (3.23) is, of course, 
developed from (3.19). Since strong kernel singularities appear when (3.23) is written for 
boundary points, once again an alternate procedure is needed to determine surface stress. 
This alternate scheme exploits the interrelationships between generalized displacement, 
traction, and stress and is the straightforward extension of the technique typically used in 
elastostatic implementation (Cruse and Van Buren, 1971). Specifically, the following can 
be obtained 

nj( = ( 3 * 24a ) 


Djfl 

2 


<r£( 0 - 


<i(0 + «ft(o) =-06ijN u ( QuZ 


^ X j N ( c\ - U N 

a/- — Q y* “»Uf 


d( “ ,JVS ' _ d( 

in which is obviously the nodal temperatures, and, 


(3.246) 

(3.24c) 


Djjki — ^^ij^ki "t 2n6ik6ji. 

Equations (3.24) form an independent set that can be solved numerically for cr£ (0 and u£)(0 
completely in terms of known nodal quantities and t^, without the need for kernel 
integration nor convolution. Notice, however, that shape function derivatives appear in 
(3.24c), thus constraining the representation of stress on the surface element to something 
less than full quadratic variation. The interior stress kernel functions, defined by (3.23), 
are also detailed in Appendix B.l. 


3.4.8 Advanced Features 

The thermoelastic formulation has been implemented as a segment of the state-of-the- 
art, general purpose boundary element computer program, GP-BEST. Consequently, many 
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additional features, beyond those detailed above, are available for the analysis of complex 
engineering problems. Perhaps, the most significant of these items, is the capability to 
analyze substructured problems. This, not only extends the analysis to bodies composed of 
several different materials, but also often provides computational efficiencies. An individual 
substructure or geometric modeling region (GMR) must contain a single material. During 
the integration process, each GMR remains a separate entity. The GMR s are then brought 
together at the assembly stage, where compatibility relationships are enforced on common 
boundaries between regions. Typically, compatibility ensures continuous displacement and 
temperature fields across an interface, however, recent enhancements to the code permit 
sliding between regions, spring contacts and interfacial thermal resistance to model air 
gaps or coating resistances. In the latter instances, discontinuities appear at the interface. 
In any case, the multi- GMR assembly process produces block-banded system matrices that 
are solved in an efficient manner. 

As another feature, a high degree of flexibility is provided for the specification of bound- 
ary conditions. In general, time-dependent values can be defined in either global or local 
coordinates. Not only can generalized displacements and tractions be specified, but also 
spring and convection boundary conditions are available. Another recent addition permits 
time-dependent ambient temperatures. A final item, worthy of note, is the availability of 
a comprehensive symmetry capability which includes provisions for both planar and cyclic 
symmetry. 

During the past two years, an interface to the well-known PATRAN graphics package 
was developed and enhanced. This interface allows the user an option to view deformed 
shapes, temperatures and stress boundary profiles or contours. A number of PATRAN- 
produced illustrations are included throughout this report. In the next section, a couple of 
examples are presented to demonstrate the validity and applicability of this boundary-only 
formulation. 
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3.5 Numerical Examples 

3.5.1 Sudden Heating of Aluminum Block 

As a first example, transient heating of an aluminum block is examined under plane 
strain conditions. The block, shown in Figure 3.1, initially rests in thermodynamic equi- 
librium at zero temperature. Then, suddenly, the face at Y = 1.0 in. is elevated to 100°F, 
while the remaining three faces are insulated and restrained against normal displacements. 
Thus, only axial deformation in the Y— direction is permitted. Naturally, as the diffusive 
process progresses, temperature builds along with the lateral stresses a xx and a zz . To com- 
plete the specification of the problem, the following standard set of material properties are 
used to characterize the aluminum: 

E = 10 x 10 6 pst, v = 0.33, 

a = 13 x 10 -6 /°F, 

k = 25in.lb./sec.in.°F, pc c = 200in.lb./in. 3o F. 

The two-dimensional boundary element idealization consists of the simple four element, 
eight node model included in Figure 3.1. A time step of 0.4 sec. is selected, corresponding 
to a non-dimensional time step of 0.5. Additionally, a finite element analysis of this same 
problem was conducted using a modified thermal version of the computer code CRISP 
(Gunn and Britto, 1984). The finite element model is also a two-dimensional plane strain 
representation, however, sixteen linear strain quadrilaterals are placed along the diffusion 
length. In the FE run, a time step of 0.2 sec. is employed. 

Temperatures, displacements, and stresses are compared in Table 3.1. Notice that the 
boundary element analysis, with only one element in the flow direction, produces a better 
time-temperature history than does a sixteen element FE analysis with a smaller time 
step. Both methods exhibit greatest error during the initial stages of the process. This is 
the result of the imposition of a sudden temperature change. Meanwhile, the comparison 
of the overall axial displacement indicates agreement to within 3% for the BE analysis 
and 5% for the FE run. A steady-state analysis via both methods produces the exact 
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answer to three digit accuracy. The last comparison, in the table, involves lateral stresses 
at an integration point in the FE model. The boundary element results are quite good 
throughout the range, however, the FE stresses exhibit considerable error, particularly 
during the initial four seconds. Actually, these finite element stress variations are not 
unexpected in light of the errors present in the temperature and displacement response. 
Recall that in the standard finite element process, stresses are computed on the basis of 
numerical differentiation of the displacements, whereas in boundary elements, the stresses 
at interior points are obtained directly from a discretized version of an exact integral 
equation. Consequently, the BE interior stress solution more nearly coincides with the 
actual response. 

3.5.2 Circular Disc 

Next, transient thermal stresses in a circular disc are investigated. The disc of radius 
initially rests at zero uniform temperature. The top and bottom surfaces are thermally 
insulated, and all boundaries are completely free of mechanical constraint. Then, suddenly, 
at time zero, the temperature of the entire outer edge (he., r — a) is elevated to unity and, 
subsequently, maintained at that level. 

The boundary element model of the disc with unit radius is shown in Figure 3.2. Only 
four quadratic elements are employed, along with quarter symmetry. Ten interior points are 
also included strictly to monitor response. In addition, the following non-dimensionalized 

material properties are arbitrarily selected for the plane stress analysis; 

E — 1.333 pc £ = 1.0 

v = 0.333 k — 1.0 

a = 0.75 

Results obtained under quasistatic conditions for a time step of 0.005 are compared, in 
Figures 3.3, 3.4 and 3.5, to the analytical solution presented in Timoshenko and Goodier 
(1970). Notice that temperatures, as well as radial and tangential stresses are accurately 
determined via the boundary element analysis. In particular from Figure 3.5, even the 
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tangential stress on the outer edge is faithfully reproduced. An extremely fine finite element 
mesh would be required to obtain a comparable level of accuracy, particularly, for the 
surface stresses. 

3.6 Summary 

A comprehensive boundary element method has been presented for transient thermoe- 
lastic analysis. This time-domain formulation requires discretization of only the surface of 
the component, and thus provides an attractive alternative to finite element analysis for 
this class of problems. In addition, steep thermal gradients, which often occur near the 
surface, can be captured more readily, since with a boundary element approach there are 
no shape functions to constrain the solution in the direction normal to the surface. For ex- 
ample, the circular disc analysis indicates the high level of accuracy that can be obtained. 
In fact, on the basis of reduced modeling effort and improved accuracy, it appears that the 
present boundary element method should be the preferred approach for general problems 
of transient thermoelasticity. 


23 



TABLE 3.1 


SUDDEN HEATING OF A CUBE 


Time 

(sec.) 

Temperature (° F ) 
at Y = 0 

Exact FE BEM 

Axial Displacement (*x in.) 
at y = 1.0 

Exact FE BEM 

Lateral Stress (ksi) 
at Y = 0.5312 
Exact FE BEM 

0.8 

4.7 

3.4 

3.8 

910 

860 

920 

- 5.6 

- 3.9 

- 5.4 

1.6 

22.0 

19.8 

20.7 

1290 

1250 

1320 

- 9.1 

- 7.7 

- 9.2 

2.4 

38.3 

36.4 

37.7 

1570 

1540 

1610 

- 11.3 

- 10.3 

- 11.7 

3.2 

51.5 

50.0 

51.5 

1780 

1760 

1840 

- 13.1 

- 12.2 

- 13.5 

4.0 

61.9 

60.7 

62.2 

1950 

1930 

2000 

- 14.4 

- 13.8 

- 14.8 

4.8 

70.1 

69.1 

70.5 

2090 

2070 

2130 

- 15.5 

- 15.0 

- 15.9 

5.6 

76.5 

75.7 

76.9 

2200 

2180 

2230 

- 16.3 

- 15.9 

- 16.7 

6.4 

81.5 

80.9 

81.9 

2280 

2270 

2310 

- 17.0 

- 16.7 

- 17.3 

7.2 

85.5 

84.9 

85.8 

2340 

2330 

2370 

- 17.5 

- 17.2 

- 17.8 

8.0 

88.6 

88.2 

88.8 

2400 

2390 

2410 

- 17.9 

- 17.7 

- 18.1 
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4. INTEGRAL FORMULATION FOR FLUIDS 

4.1 Introduction 

Attention is now shifted to the hot fluid. A number of integral formulations will 
be presented for both incompressible and compressible thermoviscous flow. In particular, 
significant effort has been directed recently toward the development and implementation of 
the convective formulations. As a result, boundary element solutions can now be obtained 
in the high Reynolds number range. 

The presentation is separated into the three classes, namely, incompressible, convective 
incompressible and convective compressible flow. Individual subsections under each head- 
ing present the governing equations, integral representations, numerical implementation 
and numerical examples. It will be evident that significant progress has been made in the 
development of boundary element methods for both incompressible cases. On the other 
hand, for the compressible case, most of the effort has been necessarily directed toward 
the derivation of new fundamental solutions, which capture the essential character of the 
flow field. 

4.2 Incompressible Thermoviscous Flow 

4.2.1 Introduction 

In the following, steady and time-dependent formulations are presented for relatively 
slow incompressible flow. The primary variables in each case are velocity, temperature, 
traction and heat flux. This is the set of variables for which boundary conditions are 
most readily defined, and for which the extension to three-dimensions is most easily ac- 
complished. As will be seen, the individual formulations have much in common. The 
major differences involve the fundamental solutions that are employed, and the treatment 
of the contributions of past events. Both formulations have been implemented within the 
computer code GP-BEST. 
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4.2.2 Governing Equations 

Application of the Principles of the Conservation of Mass, Momentum and Energy for 
an incompressible thermoviscous fluid lead to the development of the following differential 
equations: 


dvi 

dxi 


= 0 


(4.1a) 


d 2 V{ 

^ dxjdxj 


dp Dvi 
dxi ^ Dt 


4 - f% — 0 


d 2 e 

dxjdxj 



+ ^ = 0 


(4.16) 


(4.1c) 


where 

Xi Eulerian coordinate 
t time 

Vi velocity vector 
p pressure 
9 temperature 
p mass density 
/i viscosity 

k thermal conductivity 
c c specific heat 
fi body force 
tp body source, 


and the operator 


D 

Dt 


d_ 

dt 


n ‘ +V idxi 


(4.2) 


represents a material time derivative. By introducing a constant free stream velocity Ui 
and a velocity perturbation u,, such that 


Vi = Ui + Ui, 


(4.3) 
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the governing equations can be rewritten as 


du{ 

dx{ 


= 0 


(4.4a) 


d 2 u x dp dui dui dui 

<"u - Mar,- ***; h 


d 2 e 


ee .. ae de 

• dzjdxj pc ‘ Uj d^ - pc ‘ u >8^ + * = °- 


(4.46) 

(4.4c) 


Note that in equations (4.4) only the terms and pc t uj-§^ are actually nonlinear, 

although in some instances the body forces and sources may also contain nonlinearities. A 
number of distinct integral formulations are possible, depending upon which of the linear 
terms are included in the differential operator. All terms excluded from the differential 
operator, must then be grouped together as effective body forces and sources, /' and V', 
respectively. Integral formulations based upon Stokes kernels are detailed in the next 
subsection. 


4.2.3 Integral Representations 
4.2.3. 1 Steady 

In this first formulation the time-dependent terms vanish, and the entire contribution 
of the convective terms axe considered as effective body forces and sources. Thus, 


T , dui dui 


(4.5a) 


rr 86 86 
^ =- pc ^d^- pc ^8^ +1p - 


(4.56) 


As a result, the well-known fundamental solutions for incompressible Stokes flow and 
steady-state heat conduction are applicable. The integral formulation, which can be de- 
rived directly from the governing differential equation (Dargush and Banerjee, 1990c), can 
be written 

Cq/3Uq = J [Gapta — Fa0U a — G a pta]dS + J [D a pk^ka + G a pf a ]dV (4.6) 
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where 


u Q = {ui 

u 2 9} 

(4.7a) 

to = {<1 

<2 9} 

(4.76) 

fa = {fl 

h V 1 } 

(4.7c) 


axe generalized velocities, tractions, and body forces. In (4.7b), ti are the surface tractions 
defined by 

U = TijUj - pm (4.8a) 

with m representing the local unit outwaxd normal to the surface 5, and Tij the fluid 
stresses, while the heat flux is defined via 


? = 



(4.86) 


Furthermore, 


C Otp — 


Ci j 


0 

c$e 


G a p — 


Gij 0 

0 G$g 


) F a p — 


Fa 0 ] 
0 Fee 


D a pk — 


dGgp 

dx k 


a ka = [p(Uk + U k )Ui PCc(Uk + 




(4.9a, 6, c) 


(4.9d) 

(4.10a) 

(4.106) 


In the terminology of Lighthill (1952), is the momentum flux tensor or fluctuating 
Reynolds stress. Here, o° ka is labeled the generalized convective stress tensor, while t° a is 
the generalized convective traction. Both a ko and t° contain terms which axe nonlinear in 
the generalized velocities. 

In (4.9a), c;.,(£) and cge(£) are constants. When £ is inside S,Cy = and c e $ = 1. If £ is 
on the boundary then the values are determined by the relative smoothness of 5 at £. For £ 
outside the region V, both Cj j and cee are zero. Meanwhile, the kernel functions Gij,Gee, Fij 
and Fee are provided in Appendix B.2. 
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4.2. 3. 2 Time-Dependent 

For this next formulation, the effective body forces and sources are identical to those 
provided in (4.5), however, the time-dependent terms are now included in the linear oper- 
ator. The required fundamental solution for the viscous portion was first given by Oseen 
(1927), while the transient heat conduction fundamental solution is well-known (Carslaw 
and Jaeger, 1959). By applying standard methodology (Banerjee and Butterfield, 1981; 
Dargush and Banerjee, 1990d), the following governing integral equations can be derived 

Cq/JWq = J [gap *tc « — SaP *«a- 9ap * *a] dS + [d a pk * <^ka + 9aP * fa ~ 9apP^a) ( 4 - 1 D 

Note that (4.11) is similar to (4.6) for the steady case, except that Riemann convolution 
integrals over time have been introduced, along with an initial condition volume integral 
involving u° Q . Once again a° ka and t% contain terms which are nonlinear in the generalized 
velocities. Kernel functions, G Q p and F a p, developed from the instantaneous point force and 
source adjoint fundamental solutions g Q p and f Q s , are provided in Appendix B.3. It should 
be noted that these functions are considerably more complicated than the corresponding 
steady kernels. 

4.2.4 Numerical Implementation 

4.2.4. 1 Introduction 

Analytical solutions are possible for only the simplest geometries and boundary con- 
ditions. More generally, approximations must be introduced in both time and space to 
expose the practical utility of these integral equations. Consequently, in this section, state- 
of-the-art boundary element technology is applied to steady and unsteady incompressible 
thermoviscous flows. Recent boundary element developments in the fields of elastodynam- 
ics (Banerjee et al, 1986; Ahmad and Banerjee, 1988) and thermoelasticity (Dargush and 
Banerjee, 1989b, 1990a) are directly applicable for these problems. The presentation below 
will concentrate on those aspects of the numerical implementation which differ from that 
detailed in Section 3. The current implementation is limited to the two-dimensional case, 
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although certainly both of the integral formulations presented in the previous subsection 
are equally valid in three dimension. 


4. 2.4. 2 Temporal and Spatial Discretization 

For time- dependent problems, the total time interval from zero to r is subdivided into 
N equal increments of duration At. Then, the field variables i ai u a ,t%, and <T ka are assumed 


constant within each At time increment. As a result, 

N pnAr * ’ ^ 

9ap * ta — ^ / Qapdt = t Q 

n-1 n= 1 


1 

Q U Qp 


(4.12) 


with similar expressions holding for the remaining convolution integrals. This is identical 
to the treatment discussed in Section 3 for thermoelasticity. 

The methodology employed for spatial discretization of the bounding surface also fol- 
lows that described in Section 3. Thus, linear, quadratic or quartic shape functions are 
utilized to portray the functional behavior of the field variables over three-noded surface 


elements. 

However, in addition to the surface description, the domain must be discretized into 
cells in the regions where the nonlinear convective effects are important, or where nonzero 
initial conditions are present. Shape functions are once again introduced to approximate 
the geometric and functional variation with each volume cell. Thus, for any point X within 
an individual cell 

Xi{0 = M w {()x iw (4.13) 


OO = AUCKau, ( 4 - 14 ) 

where 

M Wi Mu shape functions 
Xi W nodal coordinates 
a° QU) nodal generalized convective stress . 
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The current implementation utilizes six and eight-noded cells for the geometric representa- 
tion, along with linear, quadratic, or quartic functional variation. Typical cells are depicted 
in Figure 4.1. For the quadratic cell, both serendipity (8-noded) and lagrangian (9-noded) 
variations are included. Serindipity quartic cells were found to have unsatisfactory perfor- 
mance and consequently are not available. 

As a result of the spatial discretization, the boundary integral equation for time- 
dependent thermoviscous flow can now be written 



/~*N— n+1 


N^dS — u 


n 

au> 



iriN— n+1 
F otp 


N u dS-t™ 




D^ n+l M u dV 


+ EUl/ s2pM u dV 
tel L Jv ‘ J 


n+1 


N w dS 


(4.15a) 


while for steady conditions this reduces to 

Capita — N t QuJ f GapN^dS — u QU ( F a p dS — t aw f G a pN u dS 
"i L Js m Js m Js m 

L f r 

+ E / D °0kM w dV , 

te? L Jv ' J 


(4.156) 


where M and L are the total number of surface elements and volume cells, respectively, 
and 


M 

S = ^ S m 

m= 1 


(4.16a) 


L 

v = '$2v l . (4.166) 

i=i 

The positioning of the nodal variables outside of the integrals is a key step, since now the 
integrands of (4.15) contain only known functions, which can be evaluated numerically. 

Up to this juncture, the region of interest has been assumed to be composed of a single 
volume V with surface S. However, this need not be the case. In general, space may 
be subdivided into a number of individual non-overlapping geometric modeling regions 
(GMRs). Each GMR occupies a certain volume of space, say V g , bounded by the surface 
S g . For a point £ within V g , the integration required by (4.15) need only be conducted over 
S g and V g , since the contribution to u 0 (£) from the other GMRs outside S g will be zero. 
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As a result, integration costs can be dramatically reduced by introducing multiple GMRs 
for thermoviscous flow problems. Additionally, there is no inherent requirement that all 
GMRs utilize the same physical model. For example, one GMR could employ the steady 
formulation of equation (4.6), while a second region includes the transient kernel effects 
contained in the formulation of (4.11). In any case, compatibility must, of course, be 
maintained across all GMR-to-GMR interfaces. Examples of mixed GMR formulation are 
contained in Section 4.3.6 and form the basis of the approach for fluid structure interaction 
that will be explored in Section 5. 

4. 2. 4. 3 Integration 

The evaluation of the integrals appearing in (4.15) is the next process to be examined. 
Due to the singular nature of the kernel functions G a p , F a p and D a pk considerable care must 
be exercised during numerical integration. This is particularly true for incompressible 
viscous flow, in which the final solution is extremely sensitive to errors in integration 
coefficients. In general, the integration algorithms must be much more sophisticated than 
those developed for thermoelasticity. In the present implementation, discussed in detail 
in Honkala and Dargush (1990), a number of different integration schemes are employed 
depending upon the order of the kernel singularity, the proximity of the field point £ to 
the element, and the size of the element. 

Once again consider the following three distinct categories for the surface integrals: 

(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but the kernels involve only weakly singular inte- 
grands of the In r type. 

(3) The point £ lies on the element m, and the integral has a strong £ singularity. 

In practical problems involving many elements, it is evident that most of the integration 
occurring in equation (4.15) will be of the Category (1) variety. The integrand is non- 
singular and standard Gaussian quadrature can be employed. However, for near-singular 
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causes when £ is close to element m very high order formulas axe needed to capture the 
kernel behavior. For these instances, it is beneficial to identify the point X° on the element 
nearest to and then subdivide the interval of integration about X°. Within each of 
the two subsegments a nonlinear transformation is used to further reduce the order of 
Gaussian quadrature needed for high precision. This nonlinear transformation is similar 
to that proposed by Mustoe (1984) and Telles (1987), however it should be emphasized 
that subsegmentation is still required. 

Turning next to Category (2), one finds that, unlike elasticity or potential flow, stan- 
dard Gaussian formulas alone are inadequate. Instead the terms involving In r must be 
isolated and integrated with special log-weighted Gaussian integration. The remaining 
non-singular terms comprising G a p are then evaluated utilizing standard quadrature. 

The strongly singular integrals of Category (3) exist only in the Cauchy principal 
value sense and cannot be evaluated numerically with sufficient precision. Fortunately, 
the indirect ‘rigid body’ or ‘equipotentiah method, originally developed by Cruse (1974), 
is applicable, and leads to the accurate determination of the singular block of the second 
integral in (4,15). The remainder of that integral is non-singular. Consequently, subseg- 
mentation along with standard Gaussian quadrature is adequate. 

Similar care is needed for the volume integrals, which involve the kernel D a pk con- 
taining a £-type singularity. However, for two-dimensional volume integration, this kernel 
is only weakly singular, and can be evaluated in the following direct manner. First, the 
nearest node, say A, in cell I to the point ( is determined. The cell is then subdivided 
into triangles radiating from A as shown in Figure 4.3. Next, each triangle is mapped 
onto a unit square. The apex corresponding to A is stretched to form one side of the 
square. This process essentially eliminates the ~ singularity. Finally, the square is further 
subsegmented in both radial and circumferential directions depending upon the closeness 
of f and the size of cell 1. Standard Gaussian quadrature is applied to each subsegment. 
This cell integration scheme was based on work by Mustoe (1984) for elastoplasticity. In 
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the present incompressible viscous flow implementation, tolerances have been tightened so 
that additional subsegmentation is performed, along with higher order quadrature formu- 
las. Additionally, it has been found that circumferential subsegmentation is much more 
beneficial than the radial breakup. 

In time-dependent problems, beyond the first time step, additional integration is re- 
quired. This integration involves the kernels G2p,F£ 0 and D™ 0k for n > 1. From Table 4.1, 
these are all nonsingular. As a result, a much less sophisticated integration scheme is em- 
ployed to obtain the required level of accuracy with fewer subsegments and gauss points. 
If the initial velocities are not uniform, then the nonsingular initial condition integral of 
equation (4.15a) must also be evaluated at each time step. This is accomplished in a 
manner similar to the integration of D2 0k - 

Table 4.1 - Kernel Singularities 


Kernel 


G^p for n > 1 

n, 

FqP for n > 1 

Da0k 


D 2ffk for n > 1 


Singularity Order 
In r 

non-singular 

i 

r 

non-singular 

i 


non-singular 


4. 2. 4. 4 Assembly 

Once the spatial discretization and numerical integration algorithms are completely 
defined, a system of nonlinear algebraic equations can be developed to permit an approx- 
imate solution of the thermoviscous boundary value problem. The method of collocation 
is employed by writing (4.15) at each functional mode. 
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For each time step N of a transient problem, this nodal collocation process yields 

N 

[G JV-n+1 t’ 1 - F N-r,+1 u n - G N-n+1 t or> + D N_ " + 1 <T or *] - r N u° = 0 (4-17) 

n-1 

where 

t n nodal traction vector for time step n with 3Q components 

u n nodal velocity vector for time step n with 3P components 

nodal convective traction vector for time step n with 3Q components 

a on nodal convective stress vector for time step n with 6P components 

u ° nodal initial velocity vector with 3P components 

G n unassembled matrix of size 3P x 3Q calculated from the first 

integral of (4.15) during time step n 

F n assembled matrix of size 3P x 3P calculated from the second 

integral of (4.15) during time step n, plus the c Q p contribution 
in F 1 

D n assembled matrix of size 3P x 6P calculated from the first volume 

integral of (4.15) 

pN assembled matrix of size 3P x 3P calculated from the initial condition 

integral of (4.15) 

P total number of functional nodes 

M 

Q — -4m 

m=l 

A m number of functional nodes in element m . 

All of the coefficient matrices in (4.17) contain independent blocks for each GMR in mul- 
tiregion problems. However, for any well-posed problem, the boundary conditions and 
interface relations remove all but 3P unknown components of u N and t N . Furthermore, 
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by solving (4.17) at each increment of time, all of the components of u n ,t n ,t°" and <7°" for 
n < N are known from previous time steps. Then, (4.17) can be rewritten at time NAt as 

g(x) = Ax n - DV n + G’t^ - By N 

N-l 

_ [G N-n+1 t ri - F N-r ' +1 u r> - G Ar_n+1 t°" + D N_r,+1 <7 on ] + T w u° = 0 (4.18) 


in which 

x N nodal vector of unknowns with 3P components 

y N nodal vector of knowns with 3Q components 


while A and B are the associated coefficient obtained from F 1 and G 1 . The A matrix now 
includes the compatibility relationships enforced on GMR interfaces. As a result, the GMR 
blocks in A are no longer independent, however A does remain block banded. 

The terms included in the summation of (4.18) represent the contribution of past 
events. This, along with the terms B y ;V and r A u°, can be simply evaluated once at each 
time step N with no need for iteration. Let, 

N-l 

h N _ _ By N _ y [ G N— n+l t n _ pN-n+ l u « _ Q N-n+ l,.on + D W-r>+l ff »>] + rV. (4.19) 

n=l 

Then (4.18) becomes the following nonlinear set of algebraic equations 

g(x) = Ax n - DV oW + G4 oN + = 0. (4.20) 


A closer examination of is in order. For example with N = 1 

b 1 = -By 1 +r 1 u°, (4.21a) 

while for the second time step 

b 2 = -By 2 - G 2 t' + F 2 u ! + G 2 t o1 - D 2 a o1 + r 2 u° (4.216) 

Obviously, for each step A, one new set of matrices G^.F^.D^ and T N must be determined 
via integration and assembly. Integration, particularly the volume integration needed for 
D n and T w , can be quite expensive. 
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As an alternative to the convolution approach defined above, a time marching recur- 
ring initial condition algorithm can be employed. This has been utilized by a number of 
researchers for transient problems of heat conduction, acoustics, and elasticity (Banerjee 
and Butterfield, 1981). For this latter approach, at time step N the entire contribution of 
past events is represented by an initial condition integral which utilizes u N 1 as the initial 
velocity. Thus, 

g(x) = Ax n - DV n + G 1 t oJV + b N = 0 (4.22) 

with 

b N = _ By N + r i u N-i (4.23) 

Obviously, (4.22) is identical to (4.20). Only the evaluation of b N is different. The advan- 
tage of the recurring initial condition approach is that no integration is needed beyond the 
first time step. However, volume integration is required throughout the entire domain be- 
cause of the presence of u N -‘, even for linear problems in which volume integration would 
not normally be required. 

In order to take full advantage of both methods, the present work utilizes the con- 
volution approach in linear regions, and the recurring initial condition algorithm for the 
remaining no nli near GMRs which are filled with volume cells. Since b A ' can be computed 
independently for each GMR, this new dual approach provides no particular difficulty. 

4.2.4. 5 Solution 

An iterative algorithm, along the lines of those traditionally used for BEM elastoplas- 
ticity (Banerjee and Butterfield, 1981; Banerjee et al, 1987), can be employed to solve the 
boundary value problem. However, convergence is usually achieved only at low Reynolds 
number. More generally the interior equations must be brought into the system matrix, as 
in (4.20), and a full or modified Newton- Raphson algorithm must be employed to obtain 
solutions even at moderate Reynolds number. (Similar 'variable stiffness algorithms have 
also been introduced by Banerjee and Raveendra (1987) and Henry and Banerjee (1988) 
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(4.24) 


for elastoplasticity.) Symbolically, at any iteration k, 


dg (T l \ 

S (l) 


Ax 


} = ~M 


where 


x fc+1 = x k + Ax k ( 4 -25) 

and the derivatives on the lefthand side of (4.24) are evaluated at x k . With the full 
Newton-Raphson approach, / = k and the system matrix must be formed and decomposed 
at each iteration. The out-of-core solver used in the present implementation was devel- 
oped originally for elastostatics (Banerjee et al, 1985) from the LINPACK software package 
(Dongarra et al, 1979), and operates on a submatrix level. Within each submatrix, Gaus- 
sian elimination with single pivoting reduces the block to upper triangular form. The final 
decomposed compacted form of the system matrix is stored in a direct access file for later 
reuse. Backsubstitution completes the determination of Ax fc . Iteration continues until 


ll(Ax N ) fc || 

ll(x") fc || 


(4.26) 


where € is a small tolerance, and ||x|| is the Euclidean norm of x. For the modified Newton- 
Raphson algorithm, the system matrix is not formed at every iteration, and only backsub- 
stitution is needed to determine Ax fc . 


4. 2. 4. 6 Calculation of Additional Boundary Quantities 

Once the iterative process has converged, a number of additional boundary quantities 
of interest can be easily calculated. For example, lift and drag can be calculated by numer- 
ically integrating the known nodal traction and shape function products over the surface 
elements of interest. Low order Gaussian quadrature is adequate for this integration, since 
all the functions are very well behaved. 

Furthermore, at each boundary node, the pressure?, stress <7^-, and strain rates §■£• can 
be determined by simultaneously solving the following relationships: 

VjiWnji 0 = N u «)tiu, ( 4 - 27a ) 
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(4.276) 


i(0 




dxj du{ <9A r 


5C dzj 5C 


*a( 0 


+ p(0 = o. 


(4.27c) 

(4.27rf) 


It should be emphasized that (4.27) represents a set of nine independent equations which 
are written at the boundary point £, and can be solved easily for p.cnj and at that 
point. Afterward, boundary vorticity and dilatation can be obtained, respectively, from 


Q= du 1 _d^_ 

dxi 8 x 2 

(4.28a) 

. _ dui du 2 
dx\ 8 x 2 

(4.286) 


Of course, for incompressible flow, the dilatation should be zero, but (4.28b) can be used 
as a check. 

A comprehensive PATRAN interface has also been developed. Consequently, any of 
the quantities computed above may be displayed graphically in the form of profiles or 
contours. 


4.2.5 Numerical Examples 

4.2.5. 1 Introduction 

All of the formulations discussed above have been implemented as a segment of GP- 
BEST, a general purpose boundary element code. In this section, a number of examples 
are included, primarily, to demonstrate the validity and attractiveness of the boundary 
element formulations for relatively slow incompressible flow. 

4.2. 5. 2 Converging Channel 

The two-dimensional incompressible flow through a converging channel also possesses 
a well known analytical solution which is purely radial (Millsaps and Pohlhausen, 1953). 
A comprehensive finite element study of this problem has been made by Gartling et a! 
(1977). 
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The boundary element model is shown in Figure 4.4a. The mesh contains 96 cells 
and is divided into two regions. The boundary conditions were modeled using an exact 
specification of the boundary conditions appearing in the analytical solution (Fig. 4.4a). 
Viscosity is unity, and tractions and density axe incremented to reach higher Reynolds 
numbers. The Reynolds number for this problem is defined as 

R e = PJAYAJAi (4.29) 

V 

where V 2 (Ri) is the maximum velocity in the region, which is -24.0 for the problem solved 
here. 

Figure 4.4b illustrates the results for two Reynolds numbers, indicating good accuracy 
along the entire width of the channel. Not only are the velocities accurate, but the pressures 
and tractions are very accurate also. 

It has been observed that finite element versions of this problem have several pecu- 
liarities which prevent the analytical solution from being reproduced. First of all, since 
velocities are often specified at the inlet and at the wall and centerline, ambiguous bound- 
ary condition specification results. Also, typically a parabolic “fully developed” velocity 
profile is usually specified at the inlet. However, the nonlinear solution has a flattened 
velocity distribution across the width of the channel (see Fig. 4.4b). Hence, the analyt- 
ical solution cannot be reproduced exactly if the “fully developed” profile is specified at 
the inlet. Also, the finite element modelers of this problem usually leave out the traction 
distribution at the exit and specify zero tractions there. This also gives rise to non-radial 
flow. 

The reason for so much interest in the converging flow problem is that it is one of 
the few problems possessing an analytical solution. However, by specifying a model which 
does not correspond to this problem, as in the finite element case, one cannot accurately 
compare results to the analytical solution. Any such comparisons are merely qualitative. 
In this light, the boundary element model here has utilized an exact model of the boundary 
condition and a meaningful comparison can be made. 
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4. 2. 5. 3 Transient Couette Flow 


Consider as the first transient analysis the case of developing Couette flow between 
two plates, parallel to the x-z plane, a distance h apart. Initially, both of the plates, as 
well as the fluid, are at rest. Then, beginning at time t = 0, the bottom plate is moved 
continuously with velocity V in the x-direction. Due to the no-slip condition at the fluid- 
plate interface, Couette flow begins to develop as the vorticity diffuses. Eventually, when 
steady conditions prevail, the x-component of the velocity assumes a linear profile. 

The following exact solution to this unsteady problem is provided by Schlicting (1955): 


^ n = 0 


v x (v,t) = V E er fc[2nr)i + rj\ — y; erfc[2(n + l)rji — ij\ 

n~Q 

MV’*) = 0 


(4.30a) 

(4.306) 


where 


V = 


(4 fit/py/ 2 


(4.31a, 6) 
(4.31c) 


f?1 (4/Jt/p) 1/2 

erfc(z) = 1 - erf(z) = 1 - J° e'^dy. 

All of the nonlinear terms vanish, since both v v and dv x /dx are zero. 

The two-dimensional boundary element model, utilized for this problem, is displayed 
in Figure 4.5. Four quadratic surface elements are employed, with one along each edge 
of the domain. A number of sampling points are included strictly to monitor response. 
Notice that the region of interest is arbitrarily truncated at the planes x = 0 and x = l. All 
of the boundary conditions are also shown in Figure 4.5. For the presentation of GPBEST 
results, all quantities are normalized. Thus, 


Y = | (4.32a) 

h 

T = £ (4.326) 

and the horizontal velocity is v x /V . Figure 4.6 provides the velocity profiles at four different 
times, using a time step AT = 0.025 and the convolution approach. There is some error 
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present at small times near the top plate, where the velocity is nearly zero. Results at 
Y — 0.5 versus time are shown in Figure 4.7 for several values of the time step. Obviously, 
the correlation improves with a reduction in time step and AT = 0.025 provides accurate 
velocities throughout the time history. However, even for a very large time step, the 
GPBEST solution shows no signs of instability. Error, evident in the initial portion, 
diminishes with time, and all values of AT produce the correct steady response. Further 
reduction of AT beyond 0.025 yields little benefit. Instead, mesh refinement in the y- 
direction is needed, primarily to capture the short time behavior. Figure 4.8 shows the 
GPBEST results for a model with just two, equal length, elements along each vertical side. 
The correlation with the analytical solution is now excellent. The time step selected for 
the refined model was based upon the general recommendation that 

AT = (4.33) 

c 

where is the length of the smallest element. 

The convolution approach, defined by equation (4.18), was used to obtain the results 
presented in Figures 4. 6-4. 8. Alternatively, the recurring initial condition algorithm can 
be invoked. In that case, complete volume discretization is required even for this linear 
problem. For the model of Figure 4.6, a single volume cell connecting the eight nodes is 
all that is required. The GPBEST results for different values of AT are shown in Figure 
4.9. The solutions are good for the two smaller time step magnitudes, however there is a 
slight degradation in accuracy from the convolution results. 

Interestingly, the solution in (4.30a) is identical to that for one-dimensional transient 
heat conduction in an insulated rod with one end maintained at temperature V, while the 
other remains at zero. However, in a corresponding boundary element analysis, the numer- 
ical integrations defined in (4.15a) must be calculated much more precisely for unsteady 
viscous flow than for heat conduction in order to obtain comparable levels of accuracy. 
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4.2.5. 4 Flow Between Rotating Cylinders 

As the next example, the developing flow between rotating cylinders is analyzed. The 
inner cylinder of radius rj is stationary, while the outer concentric cylinder with radius 
r a is given a tangential velocity V, beginning abruptly at time zero. The steady solution 
appears in Schlicting (1955). However, even for the transient case, the flow is purely 
circ umf erential. Thus, the governing Navier-Stokes equations reduce to 


P 


d 2 v$ 1 dv$ 

dr 2 r dr 


ve\ 

r 2 J 



dp 

dr r 


= 0 


(4.34a) 

(4.346) 


in polar coordinates (r As discussed in Batchelor (1967), separation of variables can 
be used to obtain the following solution (Honkala and Dargush, 1990) 


where 


v r (r,t) - 0 

CO 

ve(r,<) = c\r -f — -h ^ D n {J i(A n r)Yi(A n r 0 ) — Yi(A n r) J\ (A n r 0 )}e 


—A „ct 


n = l 


(4.35a) 

(4.356) 


ci = 


D„ = 


7T 


V r 0 
2 _ |*2 

2 


C 2 = —C\ r i 


An*^f ( An^t) 


T J*(\nn) — JK^r 0 ) {Yl{Xnr ° )Fln + JliXnr ° )F2n] 

Fin = — Cl [r^ J2(-^rj r o) — r i ^2(-^n r t)] + C2[./o(An r o) — •^o(^n r i)] 


F 2 n = Ci[r^V2(A n r 0 ) — r?y2(A n r<)] — C2[Y 0 (A n r 0 ) Yo(An r t)] 


(4.36a, 6) 
(4.36c) 
(4.36(f) 
(4.36c) 


and A„ is the nth root of the equation 

A(Ari)Yi(Ar 0 ) - Ji(Ar 0 )yi(Ar<) = 0. 


(4.37) 


Figure 4.10 depicts the boundary element model representing the region between the 
two cylinders. A thirty degree segment is isolated, with cyclic symmetry boundary condi- 
tions imposed along the edges 8 — 0° and 6 = 30°. The inner radius is unity, while an outer 
radius of two is assumed. Unit values are also taken for the viscosity, density and V . The 


46 



model consists of six quadratic elements and two quadratic cells. The cells, of course, are 
not needed for linear analysis utilizing the convolution approach. 

Results of the GPBEST analysis are compared to the exact solution in Figure 4.11 
for convolution and in Figure 4.12 for the recurring initial condition algorithm. In both 
diagrams, results with and without the nonlinear convective terms are plotted. The re- 
sults are quite good throughout the time history with the convolution approach, while 
some noticeable error is present at early times for the recurring initial condition solutions. 
The linear and nonlinear velocity profiles are nearly identical, as expected from the exact 
solution expressed in (4.35b). However, unlike the previous example, the nonlinear terms 
do not simply vanish from the integral equation written in cartesian form. Instead, the 
nonlinear surface and volume integrals must combine in the proper manner to produce 
the correct solution. Consequently, this problem provides a good test for the entire BEM 
formulation. 

Relative run times are shown in Table 4.2 for the different analysis types. Obviously, 
the nonlinear convolution approach is very expensive, since this involves volume integration 
at each time step. As a result, in the general implementation, convolution is only utilized 
in linear GMRs. 


Table 4.2 - Flow Between Rotating Cylinders 
(Run Time Comparisons) 


Analysis Type Time Marching Algorithm 
Linear Convolution 

Nonlinear Convolution 


Relative CPU Time 
1.0 
25.8 


Linear Recurring Initial Condition 1.5 

Nonlinear Recurring Initial Condition 1.8 
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4.2.5. 5 Driven Cavity Flow 


The two-dimensional driven cavity has become the standard test problem for incom- 
pressible computational fluid dynamics codes. In a way, this is unfortunate because of the 
ambiguities in the specification of the boundary conditions. However, numerous results 
are available for comparison purposes. 

The incompressible fluid of uniform viscosity is confined within a unit square region. 
The fluid velocities on the left, right and bottom sides are fixed at zero, while a uniform 
nonzero velocity is specified in the x-direction along the top edge. Thus, in the top corners, 
the x-velocity is not clearly defined. To alleviate this difficulty in the present analysis, the 
magnitude of this velocity component is tapered to zero at the corners. 

Results are presented for the four region, 324 cell boundary element model shown in 
Figure 4.13. Notice that a higher level of refinement is used near the edges. Spatial plots 
of the resulting velocity vectors are displayed in Figures 4.14a and b for Reynolds numbers 
(Re) of 400 and 1000, respectively. Notice that, in particular, the shift of the vortical 
center follows that described by Burggraf (1966) in his classic paper. A more quantitative 
examination of the results can be found in Figure 4.15 where the horizontal velocities on 
the vertical centerline obtained from the present GPBEST analysis are compared to those 
of Ghia et al (1982). It is assumed that the latter solutions are quite accurate since the 
authors employed a 129 by 129 finite difference grid. As is apparent, from the figure, all 
of the solutions are in excellent agreement. Finally, it should be noted that the simple 
iterative algorithm fails to converge much beyond Re = 100. Beyond that range the use of 
a Newton- Raphson type algorithm is imperative. 

In this driven cavity problem, complete volume discretization is required, since the 
nonlinear convective terms are nonzero throughout the entire domain. As a result, the 
evaluation of the volume integrals appearing in (4.6) is computationally expensive due 
to the singular nature of the kernels. Consequently, it is important to investigate the 
relative merits of a boundary element approach. To aid in this study, a finite element 
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formulation was developed based primarily on the work of Gartling et al (1977). This 
finite element implementation utilizes a penalty function approach for incompressibility, 
along with a Newton-Raphson solution algorithm. An identical sixty-four lagrangian cell 
model was selected for both the boundary element and finite element analysis. Results are 
plotted in Figure 4.16 for Re = 100. The boundary element results, though more expensive, 
are significantly more accurate. In fact, at this level of refinement, the finite element 
results show some oscillation. Clearly, for a given mesh, the boundary integral formulation 
captures more of the physics. Further comparative studies are planned for the coming 
months. 

4.2.5.6 Transient Driven Cavity Flow 

The next example involves the initiation of flow in the same square cavity. An in- 
compressible fluid of uniform density and viscosity is at rest within a unit square region. 
The velocities of the vertical sides and the bottom are fixed at zero throughout time. At 
time zero, the horizontal velocity of the top edge is suddenly raised to a value of 1000 
and maintained at that level. A gradual transition of velocities is introduced near the top 
corners to provide continuity. 

The four region, 324 cell model shown in Figure 4.13 is employed for the boundary 
element analysis. The resulting velocity vector plots at several times are shown in Figure 
4.17 for this case having a Reynolds number of 1000. The recurring condition algorithm 
was used. As in the previous two time-dependent examples, the results lead directly to 
the steady solution after a sufficient number of time steps. This steady solution correlates 
closely with the results of Ghia et al (1982), as presented in Figure 4.15. 

It should be noted that Tosaka and Kakuda (1987) have run the transient driven cavity 
at Re = 10, 000. However, their results show signs of instability even at relatively small times, 
and are compared to the steady solution of Ghia et al which also is not correct at this 
much higher Reynolds number. A valid solution in this Re range would necessitate the use 
of an extremely refined mesh, far beyond that employed by Tosaka and Kakuda or Ghia 
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et al. 


4.2.6 Summary 

The formulations presented in this section, based upon Stokes fundamental solutions, 
are suited primarily for low Reynolds number regimes. For creeping flows, all of the 
nonlinear terms vanish, resulting in a very efficient, very precise boundary-only solution. 
The resulting boundary element method is clearly superior to any of the domain based 
methods for problems of this nature, under both steady and transient conditions. 

At somewhat higher velocities, the nonlinear convective effects cannot be ignored. 
Consequently, the surface integral involving t° a and the volume integral containing a° ka in 
equations (4.6) and (4.11) are required. Since volume integration is quite computationally 
intensive, a boundary element approach becomes less attractive. This is particularly true 
when discretization is required throughout the domain, as is the case for confined flows. 
Still, for a given mesh, the boundary element formulation provides a higher degree of 
accuracy than finite difference or finite element methods, especially in the determination 
of boundary quantities. 

4.3 Convective Incompressible Thermoviscous Flow 

4.3.1 Introduction 

At high fluid velocities, the convective terms in Navier-Stokes equations tend to dom- 
inate. As a result, boundary element formulations employing Stokes kernels are inappro- 
priate, since these fundamental solutions model the effects of viscosity but not convection. 
Instead, more of the physics of the problem must be brought into the linear operator. This 
concept was clearly understood by Oseen in the early portion of the twentieth century. In 
his 1927 monograph, Oseen developed exact integral expressions for Navier-Stokes equa- 
tions using a convective fundamental solution. Unfortunately since this was well before 
the advent of the computer, he was unable to do much with his formulations beyond some 
approximate solutions at very low Reynolds number. In the present section, the work of 
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Oseen is resurrected to form the basis for an attractive boundary element method for high 
speed flows. 


4.3.2 Governing Equations 

The differential equations, governing the behavior of an incompressible thermoviscous 
fluid in the presence of a free stream velocity Ui , can be written: 

(4.38a) 


d 2 m dp dm dui , 


at 


dxjdij dii 

§* = o, 


, K 2e „ m ae , ,/ „ 
k a^-- fc ‘ Ul -^- ,lc ‘al + * =0 - 


(4.386) 

(4.38c) 


where once again represents the velocity perturbation. In (4.38), the effective body 
forces and sources axe defined as 


//= + /i (4.39a) 

ip' = —pc c Uj-^— + ip. (4.396) 

OXj 

These equations are of course identical to those presented in (4.4), except that now the 
convective terms pUjdiii/dxj and pc e UjdO/dxj are included in the linear differential operator. 
Fundamental solutions based upon (4.38) will contain the character of the flow field at 
high velocities. 


4.3.3 Fundamental Solutions 

It is instructive to begin with a look at the fundamental solution of the steady form 
of the heat equation defined above as (4.38c). In a static medium (i.e., Ui = 0), the 
fundamental solution G must satisfy 


(4 - 4o) 

in which 6 is the generalized delta function. The solution to (4.40) in two-dimensional 

space is the well-known potential flow Green’s function 

lnr 


G(x, 0 = - 


27 rk 


(4.41) 
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with 


Vi = *i~ Si 


(4.42 a) 


r - Viyi 


(4.426) 


Thus, G(x,£) represents the temperature response at x due to a unit point heat source at 
£. This response is plotted in the x\ - x 2 plane for a source at the origin in Figure 4.18. 
Radial symmetry is evident. 

However, if the medium is moving at velocity Ui, then the fundamental solution G u 
must instead satisfy 

(443) 


kg?- -*.^ + *- 0-0 

uXjdXj 0Xj 


Now, the Green’s function (e.g. Carslaw and Jaeger, 1947) is given by 

a -U k y> c / 2 k 


G u {x,() = 


27t k 


■ K ° (£)] 


(4.44) 


in which k = k/pc t . This response is plotted in Figures 4.19a-d for various magnitudes of 
an xi-directional velocity. Obviously, in a moving medium, radial symmetry is lost and 
a pronounced front-and-back effect develops. That is, at a given distance from the heat 
source, it is hottest directly downstream. 

It should be emphasized that the so-called convective fundamental solution defined in 
(4.44) actualy embodies both the processes of conduction and convection. At low velocity, 
conduction dominates producing a nearly radially symmetric response. On the other hand, 
in a high speed medium, the response is concentrated in a very narrow band downstream 
of the source. Thus, as illustrated in Figure 4.19, G v captures the transition from elliptic 
toward hyperbolic behavior. 

The corresponding convective viscous fundamental solution G ^ was first presented by 
Oseen (1911), as the solution to 

0GE 


d 2 G¥j d G% 


dxkdxk dxi 


OXk 


dGY; 


kj _ 


dx k 


0. 


(4.45a) 

(4.456) 
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The G^ tensor is given in explicit form in Appendix B.4. However, the component Gj^, 
which represents the velocity in the xj-direction due to a unit point force in the xi-direction, 
is displayed in Figures 4.20a-d. For very small Ui, the solution of (4.45) approaches the 
Stokes kernels detailed in Appendix B.2. This is shown in Figure 4.20a. Notice that, unlike 
the heat conduction response of Figure 4.19a, the static viscous fundamental solution is not 
radially symmetric. This is due to the vectorial nature of the flow, and is directly attributed 
to the t/jj/j/r 2 terms in Gij. However, as the flow velocity increases (i.e., Figures 4.20b-d), a 
stronger sense of upstream and downstream develops, and the response once again becomes 
concentrated in a narrow band ahead of the applied force. At high speed, outside of this 
band, the response is essentially zero. This behavior is not only important from a physical 
standpoint, but also can be beneficial in the development of efficient boundary element 
algorithms. 

4.3.4 Integral Representations 

The convective fundamental solutions depicted in Figures 4.19 and 4.20 capture the 
proper character of high Reynolds number incompressible thermoviscous flows, and as a 
result, can provide the basis for an attractive boundary element formulation. The corre- 
sponding integral equations, under steady conditions, can be developed directly from the 
governing differential equations (4.38). This result is, 

c o0 u a = j [G^« - Fgp- u a - G u op t u a °] dS + [Da 0 k ff ka + dV, (4.46) 

where 


a ka = [P u kUi pc t u k 0 ] (4.47a) 

<«° = **«»*• ( 4 - 476 ) 

the superscript U on the kernel functions is a reminder that these are based upon convective 
fundamental solutions. All of the kernels appearing in (4.46) are detailed in Appendix B.4. 
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In most cases the body forces, /„, are either zero or can be accounted for via a particular 
integral so that the second volume integral in (4.46) is not needed. 

In examining (4.46), it should be noted that the nonlinearities are contained in the 
surface integral involving G^/i° and the remaining volume integral, D^ 0k a^°. Specifically, 
only and <r k ° axe nonlinear, and these are both formed from the product of pertur- 
bations. For high speed flows, these perturbations are only significant in the vicinity of 
objects and in the wake. As a result, volume discretization is only needed in those areas. 
Elsewhere, the linearized Oseen approximation is adequate. 

Equation (4.46) is identical to the integral equation developed by Oseen (1927), ex- 
cept for the treatment of the nonlinear convective terms. In deriving (4.46), an additional 
integration-by-parts operation was invoked to completely eliminate the appearance of ve- 
locity gradients. 

If one is interested in the transient thermoviscous response in a medium with a more 
or less steady free stream velocity, then a time- dependent formulation is also possible. For 
this case, the time derivatives are retained in the linear operator, and the following integral 
equation results: 

C a pU 0 = J [g^0 * ta ~ fa0 * u a ~ 9ap * dS 

+ / Kp k * + 9% * U - 9 u a ppu° a \ dV ( 4 . 48 ) 

Jv 

This integral equation and the corresponding fundamental solutions have not appeared 
in the literature. The functions g v a0 axe quite involved, but can be expressed in terms of 
incomplete exponential integrals. Details will be presented next year. 

4.3.5 Numerical Implementation 

The integral representations for convective thermoviscous flow are quite similar to those 
presented in Section 4.2.3. Consequently, there is a great deal of overlap in the algorithms 
employed for their respective numerical implementation. At present, the major difference 
occurs in the schemes utilized for integration. 
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As discussed previously, the convective fundamental solutions have a much different 
character than the more familiar Stokes based kernels. The standard boundary element 
integration schemes are unable to accurately capture the localized nature of the convective 
kernels, particularly at large Reynolds number. In general, subsegmentation must be 
much more intense for singular and near-singular cases. For example, in convective near- 
singular integration, first the location X° on the element nearest to the load point ( is 
identified. Then, a graded subsegmentation pattern is defined about X° based upon criteria 
including the distance of £ to X° and the free stream velocity. For higher speed flow, 
smaller subsegments are generated. Gaussian integration order is also typically higher for 
the convective surface integration. Similar adjustments are required for volume integration 
as well. 

During this past year, some progress has been made in the development of alternate 
integration strategies for singular integration. For example, partial analytical treatment 
of the kernel has proved to be more cost effective. Also, the standard ‘rigid body’ 
technique has been extended to other known solution fields in order to indirectly calculate 
some of the singular contributions. 

However, additional effort is still needed to develop integration algorithms designed 
specifically for high speed convective kernels. In particular, the response depicted in Figure 
4.20d must be anticipated. Thus, there is no need to integrate an element which lies outside 
the narrow band of nonzero response. Furthermore, elements located partially or wholely 
within the band should be subsegmented accordingly. 

The remainder of the numerical implementation follows that discussed in Section 4.2.4. 
Thus, assembly, solution, and the calculation of additional boundary quantities are ac- 
complished in the same manner as for the Stokes kernel approach. While this is perfectly 
legitimate, full advantage has not yet been taken of the character of the convective re- 
sponse. For example, at very high speeds, as the behavior becomes hyperbolic, the system 
equations form a nearly-sequential, banded set. The present assembler and solver, which 
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were designed for elliptic systems, do not recognize this structure, and consequently, are 
quite inefficient. 

4.3.6 Numerical Examples 

4.3.6. 1 Introduction 

In order to thoroughly study the effectiveness of a boundary element approach for high 
speed flows, the above convective formulations were implemented as a segment of a state- 
of-the-art general purpose boundary element code. In the following, several numerical 
examples are presented. These examples are intended to validate the formulations, and 
to suggest the potential advantages of using a boundary element method for this class of 
problems. 

4.3. 6. 2 Burgers Flow 

The classic uniaxial linear Burgers problem provides an excellent test of the convective 
thermoviscous formulations. The incompressible fluid flows in the i-direction with uniform 
velocity U. Meanwhile, the y- component of the velocity and temperature are specified as 
U 0 and T 0 , respectively, at inlet. Both are zero at the outlet. The length of the flow field 
is L. The analytical solution (Schlicting, 1955) is 

V y = (U 0 

T = tT 0 

where 

< = {l - exp — l)]} /{l — exp[-R L }} 

with Rl = UL. 

The boundary element model employs eighteen quadratic surface elements encompass- 
ing the rectangular domain. The elements are graded, providing a very fine discretization 
near the exit, where V y and T vary substantially for large Rl • Results are shown in Figure 
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4.21 for the thermal problem and in Figure 4.22 for the viscous problem. Excellent cor- 
relation with the analytical solution is obtained in both instances for this boundary-only 
analysis, even for the highly convective case of Rl = 1000. The portion of the flow field 
just ahead of the outlet is examined more closely in Figure 4.23. The convective Oseen 
solution obviously produces a precise solution. This problem can also be solved by utilizing 
the Stokes kernels and volume cells. As seen in Figure 4.23, this latter approach is not 
quite as accurate. It should be noted that traditionally finite difference and finite element 
methods have a difficult time dealing with the convective terms present in this problem. 
Generally, ad hoc upwinding techniques must be introduced to produce stable, accurate 
solutions. On the other hand, with the convective boundary element approach the kernel 
functions contain an analytical form of upwinding. As a result, very precise BEM results 
can be obtained. 

4. 3.6. 3 Flow Over a Cylinder 

As the next convective fluids example, the oft- studied case of incompressible flow over 
a circular cylinder is considered. Initially for this problem, both the steady convective 
and non-convective formulations are utilized in the same analysis. The boundary element 
model is displayed in Figure 4.24. Note that half-symmetry is imposed.. In the inner 
region, the Stokes kernels are employed along with a complete volume discretization. Thus, 
the complete Navier-Stokes equations are represented. The outer region uses the Oseen 
kernels with a boundary-only formulation. The small non-linear contributions that would 
be present in the outer region away from the cylinder are ignored. For those more familiar 
with finite elements, each region can be thought of as a substructure or superelement. 
However, the outer region does not require a volume mesh. 

The steady-state velocity vector plot at Re = 40 is shown in Figure 4.25. The recirculat- 
ing zone, behind the cylinder, is clearly visible. Additionally, the resulting drag coefficient 
(, c D ) of 1.8 obtained from the BE analysis is within the band of experimental scatter as 
presented by Panton (1984) for the circular cylinder. 


57 


Similarly, a transient analysis can be conducted. Now a full mesh as shown in Figure 
4.26 is employed. The inner region uses a time-dependent nonlinear Stokes formulation, 
while linear Oseen kernels provide the basis for the outer infinite region. Results are shown 
in Figure 4.27a for Re = 100 at a time for which the flow is nearly fully developed. Mean- 
while, Figure 4.27b present the solution at the same time, but with a different angle of 
attack for the oncoming fluid. The results are virtually identical. This illustrates the 
relative insensitivity of boundary element solutions to the cell discretization pattern. The 
reason for this behavior, which is particularly important in modeling hyperbolic phenom- 
ena, is that so much of the boundary element formulation is analytical. Another item 
to note from these results is the completely symmetric flow patterns that were obtained. 
Asymmetry would have to be induced by perturbing either the geometry, the free stream 
velocity or the boundary conditions. 

While all of this is encouraging, the development of a simplified procedure involving 
far less volume discretization is desirable. For example, a completely linear Oseen analysis, 
which ignores all nonlinear convective terms in both regions, produces a very similar solu- 
tion, except in the vicinity of the cylinder. Vector plots from the nonlinear analysis and 
the boundary-only linear Oseen analysis are superimposed in Figure 4.28. Although it is 
difficult to distinguish between the two analyses in that plot, both produce a recirculatory 
zone behind the cylinder. Thus, the main features of the problem are captured by the 
boundary-only analysis. However, the linear solution, in general, overstates the velocities 
and velocity gradients in the neighborhood of the cylinder. Consequently, a drag coefficient 
of 3.4 is calculated, which is much higher than that found experimentally. This trend, of 
overpredicting the experimental drag, continues even to much higher Reynolds numbers 
as shown in Figure 4.29. Qualitatively, however, the behavior of the BEM Oseen solution 
is consistent with the experimental curve for Reynolds Numbers up to 100,000. 

A much improved solution can be obtained by introducing a row of cells encompassing 
the cylinder. The full nonlinear Navier-Stokes equations are solved within this inner region 
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which includes an inner and outer ring of surface elements. Exterior to the outer ring is a 
linear Oseen region. This exterior region consists simply of one matching ring of surface 
elements. Its volume extends outward to infinity, where the velocity reaches its free stream 
value. Figure 4.30 illustrates a typical mesh, along with the resulting velocity vectors. As 
Reynolds number is increased, the significant nonlinear effects concentrate nearer to the 
cylinder, so that the thickness of the inner region may be reduced. Figure 4.29 also displays 
the drag obtained by utilizing just a single row of cells. Results are quite encouraging. 

An alternative approach for high speed flows involves the conversion of the nonlinear 
volume integral into effectively a surface integral by introducing a suitable perturbation 
velocity decay function. If this is accomplished then even a nonlinear analysis reduces to a 
boundary-only solution algorithm. A concerted effort will be made in this direction during 
the coming year. 

4.3. 6.4 Flow Past Airfoils 

For illustrative purposes, a boundary-only thermoviscous analysis was conducted for 
convective flow around a pair of NACA-0018 airfoils. The boundary element model of the 
blades is shown in Figure 4.31. A hot fluid at unit temperature flows from left to right 
with a unit magnitude of the free stream velocity. Meanwhile, the airfoils are assumed to 
be stationary with their outer surface maintained at zero temperature. 

It should be emphasized that this problem was run as a boundary-only analysis, how- 
ever, a number of sampling points were included in the fluid surrounding the airfoils in 
order to graphically portray the response. First the thermal solution is examined. Figure 
4.32a depicts the temperature distribution in the fluid at a Peclet (Pe) number of ten, 
where Pe = UL/k, with fluid velocity U, thermal diffusivity »c and airfoil chord length L. 
Meanwhile, Figures 4.32b-d show the response at progressively higher Peclet number. At 
Pe = 10000, quartic surface elements were required in order to obtain an accurate solution. 
The strong convective character is quite noticeable at larger Pe as the effect of the cold 
airfoils is swept downstream. Also, in Figures 4.32c and d there is virtually no interaction 
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between the airfoils. This type of behavior is expected from a physical standpoint. It oc- 
curs in the analysis because of the banded nature of the convective fundamental solutions 
illustrated previously (e.g., Figure 4.19). However, interaction will take place if the angle 
of attack is altered. Figure 4.32e shows the response at a 30° angle for Pe = 1000. 

The velocity distribution around the airfoils follows a similar pattern. For these results 
displayed in Figure 4.33, Reynolds number is defined by Re = pU L/p. In these plots, the 
magnitude of the velocity, obtained from a boundary-only solution, is contoured. These 
results feature somewhat more interaction particularly upstream of the airfoils. It should 
be emphasized that even though a linearized solution algorithm is employed the so-called 
phenomenon of boundary layer separation can still occur. Figure 4.34 focuses on the rear 
portion of the upper blade. The contour line demarks the transition from positive to 
negative streamwise velocity, and thus very nearly identifies the point of separation. 

Next, a second row of blades is added. The modeling effort for this extension is quite 
trivial, since there is actually no discretization required beyond that needed to describe 
the airfoil surfaces. For this problem, four vertical sections of one hundred sampling points 
were included for display purposes. Velocity vectors across those sections are plotted 
in Figures 4.35 and 4.36 for Reynolds numbers of 10000 and 100000, respectively. The 
vertical spacing between the airfoils increases as one examines a through c in these 
two diagrams. The velocity profiles are noticeably affected by that spacing. However, in 
all of the plots significant velocity gradients are present. It is interesting to consider the 
level of refinement that would be necessary in a domain based finite difference or finite 
element analysis in order to capture similar gradients. 

4.3.7 Summary 

A new methodology has been presented for the solution of high Reynolds number in- 
compressible thermoviscous flow. The convective fundamental solutions that lie at the 
heart of these methods model both the diffusive character of viscosity as well as the hy- 
perbolic nature of convection. This is accomplished analytically, independent of any dis- 
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cretization pattern. Consequently, the resulting boundary element formulations are quite 
attractive particularly for higher speed unconfined flows. 

Solutions obtained for the cylinder compare favorably with experimental data. Results 
presented in Figures 4.32-4.36 for the airfoils appear to be reasonable, although these are 
not solutions to the complete Navier-Stokes equations. In particular, all terms of second 
order in the perturbation velocities have been ignored. For high speed flows, these solutions 
can be improved by including some cells in the thin boundary layer surrounding the airfoils 
and in the wake. It is not necessary to capture all of the intricacies of the flow field in 
order to obtain good engineering information on the surface of the airfoil. 

Further work is still needed in order to produce an effective analysis tool. For exam- 
ple, several promising alternatives for the representation of the nonlinear terms must be 
explored, and an intensive effort is required toward the development of efficient numerics 
tailored for the structure of the convective formulations. The latter effort should be di- 
rected toward algorithms for massively parallel machines, which provide an ideal setting 
for boundary element processing. 

4.4 Compressible Thermoviscous Flow 

4.4.1 Introduction 

Several of the previous examples have demonstrated the potential of the convective 
incompressible boundary integral formulation for flows in the high Reynolds number range. 
However, more generally, at very high speeds, compressibility of the fluid must also be 
considered. In particular, shock-related phenomena are not present in the incompressible 
formulations and kernel functions. To correct this deficiency, a compressible thermoviscous 
integral formulation is presented in this section. It should be noted that, while Oseen 
derived most of the fundamental solutions required for the incompressible case, no such 
similar solutions are available for compressibility. Consequently, considerable time and 
effort was required to derive these new approximate infinite space Green’s functions. 
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4.4.2 Governing Equations 

The conservation laws of mass, momentum and energy for a compressible thermoviscous 
fluid can be written in the following form 


dv i Dp 

+ *"° 

d 2 Vj d 2 Vi dp Dvi , 

( X + ^dxj dxi +ti dx j dx j dxi P Dt + ^~ 
, d 2 9 D6 dvi , n 

k pCf v— |-'0 = O 

dxjdxj Dt dxi 


(4.49a) 

(4.496) 

(4.49c) 


where <j> is a mass source and A is a second viscosity coefficient. All other quantities are 
defined in Section 4.2.2. Reference values for each of the primary variables are introduced 
in an effort to produce a linearized differential operator. Thus, let 


Vi — Ui + Ui 


(4.50a) 


P-Po + P 
l 6 - - 0 , o H \- 8 
P — Po + P, 


(4.506) 

(4.50c) 

(4.50(f) 


in which Ui,p o ,0 o , and p 0 are constant reference values, and Ui,p,0 and p are the perturba- 
tions. Plugging these definitions into (4.49) produces, after some manipulation, 


dui D 0 p , 

- p -eP t - dT + * =0 

d 2 uj d 2 m dp DpU j , _ 

^ +t ^dx j dx i +fi dxjdxj dxi P ° Dt + ^ 

, 9 2 e Dj „ n 

k d^dT j - p ° c ‘ Dt +tp - ° 


(4.51a) 

(4.516) 

(4.51c) 


where and ip' are now modified body mass sources, forces, and heat sources. Also, in 
(4.51), 

n ^ a 

(4.52) 


£i-l + u .± 
Dt at + 'da 
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A set of approximate fundamental solutions to (4.51) were given in the previous annual 
report (Dargush and Banerjee, 1989c). However, those solutions had two major deficien- 
cies. Firstly, the phenomenon of shock was not portrayed as expected from a physical 
standpoint, and secondly, the order of the kernel singularities was too high. 

During the current year, these fundamental solutions were abandoned. Instead, atten- 
tion was redirected toward idealizing the physical process as a combination of vortical and 
dilational motion. The vortical component is dominated by viscosity and convection, and 
is identical for both compressible and incompressible flows. On the other hand, the dilata- 
tional component must respond elastically within a convective medium. Viscous damping 
is also present. 

These considerations lead to a redevelopment of the mass conservation equation exclu- 
sively in terms of pressure. The resulting governing equations become 


c 2 d2 P +n D ° d2p D ° p + n'-o 

dxjdxj Dt dxjdxj Dt 2 

(4.52a) 

,, , N , d 7 Ui dp D 0 Ui , _ n 

0 A + p) dx,d Xl + " d Xj d Xj d Xi P ° Dt +f ' 

(4.526) 

. d 2 6 D 0 6 

t dxpir j + v =0 

(4.52c) 


where c is the speed of sound. 

4.4.3 Fundamental Solutions 

The steady two-dimensional infinite space fundamental solutions of (4.52), derived by 
Shi (1991), are presented in Appendix B.5. Since the algebraic form of these kernels is so 
complicated it is best to examine the behavior graphically. For this exercise, a forty-by- 
forty grid of sampling points was generated as shown in Figure 4.37. The source point is 
fixed at the origin, located as the central point in the grid. 

First, the component G n is plotted for various free stream velocities, expressed in 
terms of Mach number, in Figure 4.38. (Recall that Gn is the velocity in the au-direction 
at the sampling point due to a unit point force in the xi-direction at the origin.) The 
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response has some similarity to that for incompressible flow displayed in Figure 4.20. As 
the magnitude of the free stream velocity increases, a pronounced sense of flow direction 
becomes evident with the nonzero response concentrating in a narrow band ahead of the 
applied force. However, the response is always continuous, and there is a gradual evolution 
from the elliptic form at low velocity to the near-hyperbolic behavior in quickly moving 
streams. 

On the other hand, the character of G pp , representing the pressure response due to a 
unit source, is much different. At a zero Mach number, the pressure is radially symmetric 
as seen in Figure 4.39a. Increasing the Mach number to 0.9 produces a transition to 
the, by now, familiar convective form. However, at M = 1, the field suddenly becomes 
singular. Figure 4,39c shows a distinctive Mach cone at M = 1.1. It should be noted that 
the analytical kernels of Appendix B.5 produce absolutely straight lines defining the cone. 
Unfortunately, the graphics package is unable to accurately portray the discontinuity. As 
the Mach number increases further, the included angle of the cone decreases. The response 
at M = 8 is displayed in Figure 4.39d. 

Finally, Figure 4.40 shows the coupling term G ip , which measures the velocity in the 
xi-direction due to a unit source. This term also exhibits the shock-related Mach cone, 
however, now there is additionally evidence of some viscous damping of the response. 

4.4.4 Integral Representations 

The formal appearance of the governing integral equations for steady compressible 
thermoviscous flow is very similar to that provided in Section 4.3.4. Specifically, let 

C Q pU a = j [G^pta “ Fap u a] dS *f j y [G^pfa] dV 

where now 

U a = {^1 U2 P 
to = {*i t 2 dp/dn q) 

fa = if[ f2 & t//}’ 
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The major difference is, of course, in the kernel functions and F^ 0 . 

4.4.5 Summary 

New fundamental solutions were derived for compressible thermoviscous flow during 
this past year. The two-dimensional steady form is given in Appendix B.5, however solu- 
tions were also obtained for the transient case, and for three-dimensional domains. The 
contour plots of Figure 4.38 through 4.40 suggest that this latest effort has produced 
physically meaningful kernel functions. 

Although the numerical implementation of the compressible formulation has not yet 
been undertaken, a couple of characteristics of the boundary element approach should be 
noted. For high speed flows, the nonlinearities will once again be concentrated in a thin 
layer near the surface and in the wake. Thus, all of the discussion concerning high Re 
incompressible flow is valid here as well. Furthermore, with compressibility comes the 
hyperbolic phenomenon of shock. In a boundary element approach, the discontinuity can 
be captured analytically through the fundamental solution. It is not necessary to use a 
mesh to model the, generally unknown, location of the shock front. This is a distinct 
advantage for boundary elements over the domain-based methods. 
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FIGURE 4.1 TWO-DIMENSIONAL BOUNDARY ELEMENTS 
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F I.GURE 4.2 TWO-DIMENSIONAL VOLUME CELLS 
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FIGURE 4.3 - INTEGRATION SUBSEGMENTATION 
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FIGURE 4.4a - CONVERGING CHANNEL 
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FIGURE 4.4b - CONVERGING CHANNEL 
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FIGURE 4.14 - DRIVEN CAVITY FLOW 



a) Re =400 



b) Re =1000 
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HORIZONTAL VELOCITY at X=0.5 





DRIVEN CRVITY - SINGLE REGION MODEL 
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FIGURE 4.17 - TRANSIENT DRIVEN CAVITY 
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FIGURE 4.25 


FLOW OVER A CYLINDER 
VELOCITY VECTORS AT Re = 40 



FIGURE 4.26 FLOW OVER A CYLINDER - FULL MESH 
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NACA-0018 AIRFOILS • BOUNDARY ELEMENT MODEL 
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CONVECTIVE THERMOVISCOUS FLOW ( RE. PE -10. ANGLE - 0 ) 


FIGURE 4.32b 
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FIGURE 4.32e 
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CONVECTIVE THERMOVISCOUS FLOW ( RE, PE -10, ANGLE - 0 ) 


FIGURE 4.33b 
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CONVECTIVE THERMOVISCOUS FLOW ( RE, PE -1000. ANGLE - 0 ) 


FIGURE 4.33d 
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FIGURE 4.34 BOUNDARY LAYER SEPARATION 
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COMPRESSIBLE CONVECTIVE THERMOVISCOUS FLOW (M - 0.0) 


FIGURE 4.38b 
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FIGURE 4.39c 
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COMPRESSIBLE CONVECTIVE THERMOVISCOUS FLOW (M . 1.1) 


FIGURE 4.39d 
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COMPRESSIBLE CONVECTIVE THERMOVISCOUS FLOW (M . 8.0) 
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FIGURE 4.40b 
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5. FLUID-STRUCTURE INTERACTION 

5.1 Introduction 

In the previous two sections, boundary element formulations have been developed sep- 
arately for a thermoelastic structural component and for a thermoviscous fluid. However, 
the ultimate goal of this ongoing grant is to develop a single computer program to deter- 
mine the temperatures, deformation and stresses of a component exposed to a hot gas flow 
path, without the need for experimentally determined ambient fluid temperatures and film 
coefficients. While further work is still required for the fluid phase, sufficient progress has 
been made to demonstrate the utility of the overall concept. Consequently, in this section, 
problems of fluid-structure interaction will be examined. 

5.2 Formulation 

The Geometric Modeling Region (GMR) provides the vehicle for achieving interaction 
between the solid and fluid. Recall that in Section 4 different fluid formulations were 
employed in different GMRs. Now, some of the regions will use the thermoelastic solid 
boundary element model, while others utilize one of the thermoviscous fluid formulations. 
Compatibility must be enforced across all GMR interfaces, no matter which model is used 
for adjoining regions. A boundary element approach is ideal for these problems, since the 
integral equations are written directly on the interfacial surfaces. 

For demonstration purposes, consider the problem of flow past a blade as sketched in 
Figure 5.1. The blade itself is labeled GMR1, and is modeled as a thermoelastic solid. 
A boundary mesh is all that is required for this structure. Surrounding the blade is a 
thin layer of cells. This is a nonlinear thermoviscous fluid region, named GMR2, in which 
the complete Navier-Stokes equations are solved. GMR2 is enclosed by inner and outer 
surfaces composed of boundary elements. The mesh utilized for the inner surface of GMR2 
matches that employed for the blade in GMRl. Finally, the outer region GMR3, which 
extends to infinity, employs the convective Oseen kernels. The boundary element model 
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for GMR3 consists merely of the surface elements required to describe the interface to 
GMR2. Since no cells are present, the nonlinear volume and surface integrals are ignored. 
Thus, an approximation is introduced. However, as mentioned previously, outside of the 
boundary layer and wake these nonlinear contributions are negligible. (Recall that each 
region is the counterpart of a substructure or superelement commonly used in the finite 
element technology, however GMRl and GMR3 do not require any volume discretization.) 

The interface between GMR2 and GMR3 poses no particular problem. Total velocity 
and temperature from both regions are equated at each interface node, while the tractions 
and flux must be equal in magnitude but of opposite direction. The latter conditions for 
the compatibility of traction and flux are also true for the solid-fluid interface between 
GMRl and GMR2. Total temperature must, of course, be equal on this interface as well. 
However, the solid integral formulations of Section 3 are written in terms of displacement, 
while those for fluids use velocity. Consequently, a change in variable must be introduced 
to ensure complete interface compatibility. For that purpose, consider the following matrix 
form of the integral equation for a thermoviscous fluid: 

i4 0 lT 


0 coe 


{ v * \ _ 0 f U 1 _ Fij 0 / \ . { Rj X (5 i\ 

\0J [o Ge$\ 1 q J [0 F ee \ \ 0 j ^ \R e ) ’ [ } 


The contributions from nonlinearities and past time steps are all contained in Rp, as are 
any terms associated with the translation from perturbed velocity to total velocity d*. 
Meanwhile, a similar expression written for a thermoelastic solid becomes 


C{j 0 

0 coe 


f Ui 1 _ Gij 0 f U 1 _ Fij 0 / Ui \ i / X /c o\ 

10 J L g *3 G es\ \qj [ F ej F ee { 6 } ^ \ Re j ' K } 


where Ui is the total displacement. This must be rewritten in terms of total velocity t 
where 

dui 


Vi = 


dr ' 


( 5 . 3 ) 


After invoking properties of the convolution integrals that are present in the original inte- 
gral equation (3.2), the appropriate representation for the solid can be written 


Cii 


0 

0 coo 
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Gij 
[ Gej Gee J 
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in which Gij,G e j and F$j are now modified kernel functions and Rp is the corresponding 
right-hand-side contribution. However, at this point, the fluid formulation (5.1) and the 
solid formulation (5.4) are completely compatible, and are in an ideal form to solve quite 
general interaction problems. 

5.3 Numerical Implementation 

The boundary element code, GPBEST, was generalized so that any combination of 
solid and fluid regions could be accommodated. Also, the modified thermoelastic kernels 
of equation (5.4) were implemented. The entire GPBEST input is free format and keyword 
driven. Output is provided on a region-by-region basis, and thus contains only informa- 
tion pertinent to the region type. Displacements, temperatures, stresses and strains are 
detailed for solid GMRs, while velocities, temperatures, stresses, pressures, strain rates 
and vorticities are output for fluid regions. In all cases, a complete PATRAN interface is 
available, so that any quantities can be plotted. 

5.4 Numerical Examples 

5.4.1 Introduction 

In this subsection a couple of examples will be presented to highlight the attractiveness 
of the present coupled boundary element approach. Flow past a thick-walled cylinder and 
an airfoil are considered. Both steady and transient conditions are examined, and a number 
of additional features of the GP-BEST implementation are explored. 

5.4.2 Steady Response of a Thick Cylinder 

For the first example, a thick-walled stainless steel cylinder rests under plane strain 
conditions in a stream of hot gas. The cylinder has an outer diameter of 1.0 in. and a 
thickness of 0.125 in. The inner surface of the cylinder is maintained at a temperature of 
0°F, while the gas temperature in the free stream is 1000°F. The following thermoelastic 
properties are assumed for the solid cylinder 
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F = 29. x 10 6 psi, v - 0.30 

a ~ 9.6 x 10 _6 in./in°F 
k — 6.48 in. lb. /sec. in °F 

p = 7.34 x 10" 4 lb.sec. 2 /in 4 c £ = 3.83 x 10 5 in.lb.in./lb.sec. 2o F. 

Additionally, the thermoviscous properties of the hot gas are taken as 
p = 5.30 x 10“ 9 lb. sec/in. 2 

k = 7.28 x 10~ 3 in.lb./sec.in°F 

p = 3.69 x 10" 8 lb.sec. 2 /in 4 c p = 9.49 x 10 5 in.lb.in./lb.sec. 2o F. 

Fluid velocities of 144 in./sec., 1440 in./sec. and 14400 in./sec., corresponding to Reynolds 
Numbers of 10 3 , 10 4 and 10 5 , are examined. In all cases, the hot gas flows from left to right, 
and only the steady response is considered. 

At Re = 1000, the maximum temperature in the cylinder is only 98°F, and the peak 
compressive axial stress is 36 ksi. However, when the fluid velocity is increased to attain 
an Re = 10, 000 a much more significant response is obtained. The temperature contours 
are shown in Figure 5.2a, the deformed shape is depicted in Figure 5.2b, and Figure 
5.2c illustrates the axial stress distribution. It should be noted that in Figure 5.2b the 
deformation has been scaled by a factor of 100. The effects of convection are quite evident 
in all three diagrams. With Reynolds number increased to 100,000 these effects become 
even more pronounced, as seen in Figures 5.3. Now the peak metal temperature has 
reached 918°F. 

5.4.3 Airfoil Exposed to Hot Gas Flowpath 

In this final example, an NACA0018 airfoil with an internal cooling passage is exposed 
to the flow of a hot gas. The boundary element model for the airfoil is shown in Fig- 
ure 5.4. Each dash represents an individual quadratic surface element. Throughout this 
problem, the outer gaseous region is modeled as a linear steady convective domain. Thus, 
a boundary-only exterior GMR is employed for the fluid. The hot gas at 1000°F flows 
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from left to right, while the inner surface of the airfoil is maintained at 200°F. Material 
properties from the previous example are once again used to characterize both the solid 
and fluid. 

For the first set of investigations, the behavior of the airfoil is determined under steady- 
state conditions. Figure 5.5a displays the deformed shape at a Reynolds number of 1000 
(based upon chord length). The solid line represents the final deformed shape, except 
that displacements have been scaled by a factor of twenty-five. Meanwhile, Figures 5.5b 
and c present the profiles of temperature and axial stress, respectively, along the upper 
surface of the airfoil. At this relatively slow speed flow, the airfoil is only effected near 
its leading edge. More significant response is shown in Figures 5.6a-c for Re = 10,000 and 
Figures 5.7a-c for Re = 100,000. In the latter case, the temperature at the stagnation point 
is nearly that of the free stream. All three cases considered so far have assumed an angle 
of attack of 0° with respect to the x-axis. Consequently, the response of the upper and 
lower surfaces is identical. Next, the angle of attack (a) is modified to 5° and 10°. Results 
for these cases are shown in Figures 5.8 and 5.9, respectively. Considerable asymmetry 
between upper and lower surfaces is now evident, although peak values of temperature and 
stress are essentially unaffected. 

Thermal barrier coatings are often employed to reduce the metal temperatures and 
stresses in hot section components. The benefit of such coatings can easily be evaluated 
with the present boundary element formulation. Consider, for example, a coating material 
with thermal conductivity k = 0.50 in.lb./sec.in.°F sprayed to a thickness of .0095in. This is 
equivalent to an interfacial thermal resistance of .021 sec.in°F/in.lb., which can be specified 
on the fluid-to-solid GMR interface. Results are displayed in Figure 5.10 for Re = 100,000 
at a = 10°. Peak airfoil temperature is reduced from 976°F to 738°F by introducing this 
particular thermal barrier coating. 

Finally, it is of considerable interest to examine the transient response of the airfoil. 
At time zero, the airfoil is in thermal equilibrium at a temperature of 200°F. Suddenly, 
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it is subjected to the hot gas stream with Re = 100,000 and a = 10°. The response of the 
upper surface at 1 msec., 2msec., 5 msec., and 10 msec, is shown in Figures 5.11-5.14. 
For this transient case, the peak stress occurs slightly offset from the tip of the airfoil. 
Additionally, the stress a yy reaches a maximum at approximately 2 msec., while a zz and 
the temperature continue to climb to their steady-state values. This is true of the axial 
stress only because of the assumption of plane strain. In a full three-dimensional analysis, 
g zz would also have a higher peak during transient state. 
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FIGURE 5.1 - FLUID-STRUCTURE INTERACTION CONCEPTUAL MODEL 
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AIRFOIL - BOUNDARY ELEMENT MODEL 
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FIGURE 5.6 - AIRFOIL (STEADY; Re = 10,000; a = 0°) 
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FIGURE 5.9b-e - AIRFOIL (STEADY: 
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6. BEM FOR RELATED PHYSICAL PHENOMENA 

During the course of the investigation of the hot fluid- structure problem, a number of 
related technologies have been opened to analysis by the boundary element method. In this 
section, several of these potential applications are discussed. Most of the advancements 
depend upon the development of new fundamental solutions. For each case, a systematic 
procedure can be applied to obtain the required fundamental solution. This same procedure 
was developed and refined during the derivation of all of the kernel functions presented m 

Sections 3 and 4. 

Perhaps the most interesting of these applications involve either moving sources or 
moving media. An example of the former kind is the determination of residual stresses m 
welds. As part of the NASA/HOST program, the boundary element code BEST3D was 
developed for the inelastic analysis of structures. Included m that code are a number of 
elastoplastic and viscoplastic material models that would be suitable for the weld problem. 
However, the temperature in the weld and adjoining structure is not known a priori, and 
a transient heat conduction analysis is required which accounts for the speed of the weld. 
The desired integral formulation for this thermal analysis is quite similar to that discussed 
for convective flow in Section 4. In addition, the fundamental solution that is needed for 
moving heat sources has already been derived as part of the present work. The other 
major advancement in boundary element technology that is required to solve the weld 
problem involves the development of more sophisticated nonlinear solution algorithms. It is 
envisioned that the modified Newton-Raphson schemes, employed for thermoviscous fluids, 
will provide the basis for that development. It should be noted that similar problems, such 
as frictional heating, grinding, and machining could also be studied utilizing the moving 
heat source approach. 

The hot viscous fluid formulations presented in Section 4 are quite general, and conse- 
quently, applicable to a wide range of physical processes. For example, the incompressible 
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integral equations could be used to solve the flow problem in injection molds, or the con- 
vective formulations could be applied to investigate the cooling of electronic components. 
Furthermore, some relatively minor extensions would provide significant benefits. The in- 
clusion of a buoyancy term based upon the Boussinesq approximation, would permit the 
examination of the thermally-induced flow in lakes or the slow heating of a room. The 
addition of an extra equation involving the concentration of a diffusing substance provides 
the opportunity to investigate the spread of pollutants in a convective environment. 

As mentioned previously, once the techniques for obtaining fundamental solutions have 
been mastered, a wide range of physical phenomena can be analyzed via boundary element 
approach. Recent work by Kaynia and Banerjee (1990) has focused on the development 
of fundamental solutions for dynamic poroelasticity. These solutions will be utilized in 
a BEM (Chen, 1991) for the analysis of soil-structure interaction under seismic loading. 
The analogous problem of dynamic thermoelasticity, which includes the important case of 
thermal shock, can also be solved with the same formulation. 

The coupling approach discussed in Section 5 can be used not only to solve the ther- 
moviscous fluid-structure problem, but also to investigate flutter. In this case, frequency- 
dependent formulation solutions are required. The infinite space solution for periodic 
elastodynamics of solids is well-known (Banerjee and Butterfield, 1981), while that for a 
linearized Oseen fluid could be derived. The frequency domain BEM analysis would be an 
extension of the work done for the NASA/HOST program and contained in BEST3D. 

There currently exists no satisfactory numerical nor analytical techniques to effectively 
deal with all of the physical phenomena mentioned in the preceding paragraphs. However, 
as an indirect result of the present hot fluid-structure grant, boundary element formulations 
and implementations are now possible for each case. 
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7. SUMMARY 

A new methodology has been developed for hot fluid structure interaction based upon 
an integrated boundary element approach. As a part of this effort, significant advances 
have been in the analysis of both the solid and the surrounding fluid. 

Section 3 detailed a boundary-only, time domain formulation for the analysis of ther- 
moelastic solids. Not only does this approach eliminate the need for volume discretization, 
it also permits the accurate determination of surface temperatures and stresses which are of 
primary interest in hot section components. Thus, this boundary element method is a suit- 
able substitute for finite elements for this entire class of problems. The two-dimensional 
formulation was presented here, however three-dimensional and axisymmetric methods 
have also been developed. 

As mentioned previously, most of the effort expended during this research program 
has been directed toward the development of appropriate boundary element methods for 
thermoviscous fluids. This was necessary because only rudimentary formulations were 
available in the literature. For slow creeping flows a boundary-only method was developed 
for both steady-state and transient problems. In these flows, the nonlinear convective 
terms are negligible. As the fluid velocity is increased to moderate levels, these convective 
effects can no longer be ignored. Consequently, volume discretization is required and 
the boundary element approach based upon Stokes fundamental solutions becomes less 
attractive primarily due to the cost of cell integration. However, it should be noted that 
the resulting boundary element solutions are typically very accurate. 

At higher speeds, when the convective effects dominate the entire problem, it no longer 
makes sense to use the viscous-based Stokes kernels. Instead, Oseen convective fundamen- 
tal solutions are employed. As demonstrated in Section 4.3, those new kernels embod} 
more of the physics of high Reynolds number flows. In fact much of the character of the 
problem can be captured with a linear boundary-only analysis. However, if more accuracy 
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is desired, volume cells can be added to the nonlinear portions of the flow field; namely, 
the thin boundary layer and the wake. These cells basically are used only to correct the 
linear solution. It is generally not necessary to capture all of the minute details of the 
flow in order to obtain the desired surface information, although, for example, there is no 
reason that turbulence models could not be introduced within the nonlinear regions. 

For compressible flows, the corresponding fundamental solutions do not appear in 
the published literature, and, consequently, had to be developed. A new set of kernel 
functions, derived during this past year, were presented in Section 4.4. As shown in the 
diagrams, these kernels explicitly contain the hyperbolic nature of shock. Consequently, 
the boundary element formulations, based upon these Green’s functions, will be able to 
model the shock front without the need for a discretization pattern. This will provide a 
significant advantage for the boundary element approach over any of the existing numerical 
techniques. 

Finally, in Section 5, the boundary formulations for a thermoelastic solid were com- 
bined with those of a thermoviscous fluid to create a novel hot fluid structure interaction 
capability. Since integral equations are written directly on the fluid-structure interface, 
the BEM approach is ideally suited for this class of problems. A couple of examples were 
included to demonstrate the attractiveness of this method in terms of model generation 
and results interpretation. Additionally, it should be emphasized that all of the numerical 
solutions included in this report were obtained on a standard desktop SUN SPARCstation 
1 . 

In light of all of the above developments, it must be concluded that an effective new 
approach has been identified for computational fluid dynamics and hot fluid-structure 
interaction. However, much additional effort is needed. Some of the required tasks are 
outlined in the next section, which defines the future direction of our research effort. 
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8. FUTURE DIRECTION 

Despite the progress that has been made during the course of this research program, 
the present boundary element approach for hot fluid-structure interaction is still primar- 
ily limited by the ability to properly model and efficiently calculate the response of the 
surrounding fluid. The boundary element formulations for fluids are particularly attrac- 
tive at the two extremes of low and high speed flows. At low velocities, the extensive 
boundary element technology developed for solids is directly applicable since the problems 
are primarily elliptic. In the intermediate range, it is quite appropriate to consider the 
combination of methods, with finite element or finite difference methods employed in the 
nonlinear regions and boundary elements for the outer linear portions of the flow field. 
Some attention will be given to this approach in the coming year. 

However, for high speed flows, the character of the response changes. Nonlinearities 
axe confined to the vicinity of the structure and the behavior becomes more hyperbolic. 
Consequently, a purely boundary element approach once again becomes most attractive. 
However, the necessary integration, assembly, and solver technologies have never been de- 
veloped for this type of system. Instead of the standard family of volume cells, decay func- 
tions should be introduced to effectively reduce the volume integration to a surface-based 
computation. Furthermore, during all integration, the banded nature of the fundamental 
solutions should be recognized. Similarly, efficiencies can be introduced during assembly, 
where currently many zeroes are processed. In the solver, advantage must be taken of the 
nearly- sequential structure of the system equations. The implementation of these ideas 
would result in a very efficient method for high Reynolds number flows, particularly in a 
massively parallel computing environment. In fact, the fluid dynamics boundary element 
algorithm, including the features outlined above, is ideally suited for that environment, 
since it involves a large number of computationally-intensive independent processes. 

Additionally, further work is needed on the implementation of the convective compress- 
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ible fundamental solutions, and the corresponding three-dimensional formulations must be 
developed. A number of planned research activities for the coming years are listed below, 
primarily in chronological order. Of course, the amount of progress that can be achieved 
in 1991 will be largely dependent on the level of funding. 

Research Plan 

• Development of a nonlinear boundary layer representation in terms of decay functions. 

• Implementation of the new convective compressible formulation. 

• Development of revamped convective integration algorithms based upon the nature of 
the kernels. 

• Investigation of several more realistic problems of hot fluid-structure interaction. 

• Incorporation of a nonlinear finite element region. 

• Restructuring of the assembly routines for convective flows. 

• Development and implementation of three dimensional formulations for fluid-structure 
interaction. 

• Development of a nearly-sequential banded solver. 

• Development of massively parallel computing algorithms. 


143 


APPENDIX A - References 

Ahmad, S. and Banerjee, P.K. (1988), ‘Transient Elastodynamic Analysis of Three Dimen- 
sional Problems by BEM,’ Int. Jour. Numerical Methods in Engineering, Vol. 26, No. 8, 
pp. 1560-1580. 

Baneriee P.K., Ahmad, S. and Manolis, G.D. (1986), ‘Transient Elastodynamic Analysis 
of Three-dimensional Problems by Boundary Element Method, Earthquake Engineering 
and Structural Dynamics, Vol. 14, pp. 933-949. 

Banerjee, P.K. and Butterfield, R. (1981), ‘Boundary Element Methods in Engineering 
Science,’ McGraw-Hill, London. 

Banerjee, P.K. and Morino, L. (1990), Boundary Element Methods in Nonlinear Fluid 
Dynamics, Developments in Boundary Element Methods-6, Elsevier Applied Science, Eng 
land. 

Banerjee, P.K. and Raveendra, S.T. (1987), ‘A New Boundary Element Formulation for 
Two-dimensional Elastoplastic Analysis,’ Jour, of Engrg. Mech., ASCE, V113, No. 2, pp. 
252-265. 

Banerjee, P.K., Wilson, R.B. and Miller, N (1985k ‘Development of a Large BEM Sys- 
tem for Three-dimensional Inelastic Analysis ’ in Ad^^ Topi^ 

Analysis, ed. T.A. Cruse, A.B. Pifko and H. Armen, AMD-V72, ASME, New York. 

Banerjee, P.K., Wilson, R.B. and Miller, N. (1988), ‘Advanced Elastic and Inelastic Three- 
dimensional Analysis of Gas Turbine Engine Structures by BEM, Int. J. Num. Meth. 
Engrg., V26, pp. 393-411. 

Banerjee, P.K., Wilson, R.B. and Raveendra, S.T. (1987), ‘Advanced Applications of BEM 
to Three-dimensional Problems of Monotonic and Cyclic Plasticity, Int. Jour. Mech. 
Sciences, V29, No. 9, pp. 637-653. 

Batchelor, G.K. (1967), An Introduction to Fluid Dynamics, Cambridge University Press, 
Cambridge, U.K. 

Boley, B.A. and Weiner, J.H. (1960), ‘Theory of Thermal Stresses,’ John Wiley and Sons, 
New York. 


Burggraf O.R. (1966), ‘Analytical and Numerical Studies of the Structure of Steady Sep- 
arated Flows,’ J. Fluid Mech., V24, Part 1, pp. 113-151. 

Carslaw, H.S. and Jaeger, J.C. (1959), Conduction of Heat in Solids, Clarendon Press, 
Oxford. 

Chaudouet, A. (1987), ‘Three-dimensional Transient Thermoelastic Analysis by the BIE 
Method,’ Int. J. Num. Meth. Engrg., V24, pp. 25-45. 

Chen, J. (1991), Dynamic Poroelasticity Via the Boundary Element Method, forthcoming 
Ph.D! Dissertation, State University of New York at Buffalo. 

Cruse T A 11974') ‘An Improved Boundary Integral Equation Method for Three Dimen- 
S’Elastic Stress Analysis,’ Comp, and Struct., V4, pp. 741-754. 


144 


Cruse, T.A. and VanBuren, W. (1971), ‘Three-dimensional Elastic Stress Analysis of a 
Fracture Specimen with an Edge Crack/ Int. J. Fract. Mech., V7, pp. 

Cruse, T.A., Snow, D.W. and Wilson, R.B. (1977), ‘Numerical Solutions in Axisymmetric 
Elasticity,’ Comp, and Struct., V7, pp. 445-451. 

Dargush G.F. (1987), BEM for the Analogous Problems of Thermoelasticity and Soil 
Consolidation, Ph.D. Dissertation, State University of New York at Buffalo. 

Dargush, G.F. and Banerjee, P.K. (1988), ‘Development of an Integrated BEM for Hot 
Fluid-Structure Interaction,’ Advanced Earth-to-Orbit Propulsion Technology Conference, 
NASA CP-3012, Huntsville, May 1988. 

Daxgush, G.F. and Banerjee, P.K. (1989a), ‘Development of an Integrated BEM for Hot 
Fluid-Structure Interaction,’ International Gas Turbine and Aeroengine Congres 
position, ASME, Paper 89-GT-128, Toronto; also m J. Eng. Gas Turbines and Pover, 

V112, pp. 243-250. 

Daxgush, G.F. and Banerjee, P.K. (1989b), ‘Development of a B °und^y Element Method 
for Time-dependent Planar Thermoelasticity, Int. J. Solids Struct., V25, pp. 999 1021. 

Daxgush, G.F. and Banerjee, P.K. (1989c), Development of an Integrated BEM Approach 
for Hot Fluid Structure Interaction, NASA Annual Report, Grant NAb3-tl2. 

Dargush, G.F. and Banerjee, P.K. (1990a), ‘Boundary Element Methods in Three Dimen- 
sional Thermoelasticity,’ Int. J. Solids Struct., V26, pp. 199-216. 

Daxgush, G.F. and Banerjee, P.K. (1990b), ‘Time Dependent Axisymmetric Thermoelastic 
Boundary Element Analysis/ submitted to Int. J. Num. Meth. Eng. 

Daxgush, G.F. and Banerjee, P.K. (1990c), ‘Advanced Boundaxy Element Methods for 
Steady Incompressible Thermoviscous Flow,’ in Developments m BEM-6, ed. P.K. Baner- 
jee and L. Morino, Elsevier Applied Science Publishers. 

Dareush G.F. and Banerjee, P.K. (1990d), ‘A Time-dependent Incompressible Viscous 
BEM for Moderate Reynolds Number,’ in Developments in BEM-6, ed. P.K. Banerjee an 
L. Morino, Elsevier Applied Science Publishers. 

Daxgush, G.F. and Banerjee, P.K. (1991a), ‘A Boundaxy Element Method for Steady 
Incompressible Thermoviscous Flow/ to appear in Int. J. Num. Meth. Eng. 

Daxgush, G.F. and Banerjee, P.K. (1991b), ‘A Time Dependent Incompressible Viscous 
BEM for Moderate Reynolds Numbers/ to appear in Int. J. Num. Meth. Eng. 

Daxgush, G.F., Banerjee, P.K. and Dunn, M.G. (19871, Development of an ^grated BEM 
Approach for Hot Fluid Structure Interaction, NASA Annual Report, Grant NAG3-712. 

Daxgush, G.F., Banerjee, P.K. and Honkala, K.A. (1988), Development of ^ Integrated 
BEM Approach for Hot Fluid Structure Interaction, NASA Annual Report, Grant NAG3 

712. 

Deb, A. and Banerjee, P.K. (1989), ‘A Comparison Between Isoparametric Lagrangian 
Elements in 2D BEM,’ Int. J. Num. Meth. Eng., V28, pp. 1539-1555. 


145 


Dongarra, J.J. et al (1979), Linpak User’s Guide, SIAM, Philadelphia. 

Gartling, D.K., Nickell, R.E., Tanner, R.E. (1977), ‘A Finite Element Convergence Study 
for Accelerating Flow Problems,’ Int. J. Num. Methods Eng., Vll, pp. 1155-1174. 

Ghia, U., Ghia, K.N. and Shin, C.T. (1982), ‘High-Re Solutions for Incompressible Flow 
Using the Navier-Stokes Equations and a Multignd Method, J. Comp. Physics, V48, pp. 
387-411. 

Gladden, H.J. (1989), ‘Aerothermal Loads on Actively Cooled Components: Analyses and 
Experiment,’ HITEMP Review, NASA Conference Publication 10039, Cleveland, Ohio, 
Oct. 31-Nov. 2, pp. 68.1-68.12. 

Gunn, M.J. and Britto, A.M. (1984), CRISP User’s and Programmer’s Guide, Engineering 
Department, Cambridge University. 

Henry, D.P. and Banerjee, P.K. (1988), ‘A Variable Stiffness Type Boundary Element 
Formulation for Axisvmmetric Elastoplastic Media,’ Int. Jour, for Num. Methods 
Engrg., V25, pp. 1005-1027. 

Honkala, K.A. and Dargush, G.F. (1990), ‘Transient Problems of Incompressible Thermo- 
viscous Flow via BEM,’ submitted to Jour. Comp. Physics. 

Kaynia, A.M. and Banerjee, P.K. (1990), ‘Fundamental Solutions of Biot’s Equation of 
Dynamic Poroelasticity,’ submitted to J. Appl. Mech. 

Latchat, J.C. and Watson, J.O. (1976), ‘Effective Numerical Treatment of Boundary Inte- 
gral Equations: A Formulation for Three-dimensional Elastostatics, Int. J. Num. Metii. 
Engrg., V10, pp. 991-1005. 

Millsaps, K. and Pohlhausen, K. (1953), ‘Thermal Distributions in Jeffery-Hamel Flows 
Between Nonparallel Plane Walls,’ Journal of the Aeronautical Sciences, March, pp. 18/ 

196. 

Mustoe G.G.W. (1984), ‘Advanced Integration Schemes Over Boundary Elements and 
Volume Cells for Two- and Three-dimensional Nonlinear Analysis, in Developmen s m 
Boundary Element Methods - III, ed. P.K. Banerjee and S. Mukherjee, Applied Science 
Publishers, England. 


Oseen, C.W. (1911), Uber die Stokes’sche Formel und iiber eine verwandte Aufgabe in der 
Hydrodynamik II, Ark. f. mat., astr. och fysik, V7. 

Oseen, C.W. (1927), Neuere Methoden und Ergebnisse in der Hydrodynamik, Akad. Ver- 
lagsgellschaft , Leipzig. 

Panton, R.L. (1984), Incompressible Flow, John Wiley and Sons, New York. 

Piva, R. and Morino, L. (1987), ‘Vector Green’s Function Method for Unsteady Navier- 
Stokes Equations,’ Meccanica, Vol. 22, pp. 76-85. 

Piva, R. Graziani, G. and Morino, L. (1987) ‘Boundary Integral Equation Method for Un- 
steady Viscous and Inviscid Flows,’ IUTAM Symposium on Advanced Boundary Element 
Method, San Antonio, Texas. 


146 



Rizzo, F.J. and Shippy, D.J. (1977), ‘An Advanced Boundary Integral Equation Method 
for Three-dimensional Thermoelasticity, ’ Int. J. Num. Meth. Eng. Vll, pp. 1753-1768. 


Schlicting, H. (1955), Boundary Layer Theory, McGraw-Hill, New York. 

Sharp, S. and Crouch, S.L. (1986), ‘Boundary Integral Methods for Thermoelasticity Prob- 
lems,’ J. Appl. Mech., V53, pp. 298-302. 

Shi Y. (1991), Boundary Element Methods for Compressible and Incompressible Viscous 
Flow, forthcoming Ph.D. Dissertation, State University of New York at Buffalo. 

Stroud, A.H. and Secrest, D. (1966), Gaussian Quadrature Formulas, Prentice Hall, New 
York. 

Telles J.C.F. (1987), ‘A Self-Adaptive Co-ordinate Transformation for Efficient Numerical 
Evaluation of General Boundary Element Integrals,’ Int. J. Num. Meth. Engrg., V24, pp. 
959-973. 

Timoshenko, S.P. and Goodier, J.N. (1970), Theory of Elasticity, McGraw-Hill, New York. 

Tosaka, N. and Kakuda, K. (1986), ‘Numerical Solutions of Steady Incompressible Viscous 
Flow Problems by Integral Equation Method,’ pp. 211-222 in.R.P. Shaw et al, eds. Proc. 
4th Inti. Symp. Innov. Num. Methods Engrg., Springer, Berlin. 

Tosaka, N. and Kakuda, K. (1987), ‘Numerical Simulations of Laminar and Turbulent 
Flows by Using an Integral Equation,’ Boundary Element IX, eds. Brebbia, Wendland 
and Kuhn, pp. 489-502. 

Tosaka, N. and Onishi, K. (1986), ‘Boundary Integral Equation Formulations for Unsteady 
Incompressible Viscous Fluid Flow by Time- differencing, Engineering Analysis, , o. 
2, pp. 101-104. 


147 


APPENDIX B.l - 
Kernels for Thermoelasticity 


This appendix contains the detailed presentations of all the kernel functions utilized in 
the formulations contained in Section 3. Two-dimensional (plane strain) kernels are pro- 
vided, based upon continuous source and force fundamental solutions. For time-dependent 
uncoupled quasistatic thermoelasticity the following relationships must be used to deter- 
mine the proper form of the functions required in the boundary element discretization. 

That is, , , 

G2 p (X-Q = G 0 p{X-t,n At) for n = 1 

GSp(* - 0 = G aP (X - (, nAt) - G a pX - £, (n - 1)A<) for n > 1, 
with similar expressions holding for all the remaining kernels. In the specification of these 

kernels below, the arguments (X — £,t) are assumed. The indices 

i,j,k,l vary from 1 to d 

a,/3 vary from 1 to (d + 1) 

$ equals d + 1 

where d is the dimensionality of the problem. Additionally, 
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& coordinates of field point 
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For the displacement kernel, 

Gio — 0 



whereas, for the traction kernel, 


^ tJ 47rr(l 


J ^ViyjVknk ^ _ ^ijVkrik + yjTij ^ ^ 


4- 




F ie = 0 

5 ( 1 * 5 ) [(*5?) am- mmA 


In the above, 


T? = Jctyfi 

~ _L 

” pc £ 


•M-r 
iM = ? ( 



1 - e-^ 4 ) 


ff4(>?) = 


§s(»?) = 




/e(»?) = 


M*?) = 


2 2 


/<(»*) = 


For the interior stress kernels, 

r. _ 2jxz/ c dGpi t ( dGpi &Gpj\ oc 

Epij - + M VW ~WJ~ P ii pe 

n 2^1/ f aF/ji t /'SFei dFpj \ at 

Dpij - — to 6ij ^T + /i V"^7 + ^rj“ ^ ijF/3S 


where 


aGv 


a^fc 8flr fi(i — u) 




VkVk _ 6jkVi _ 6ikVj \ 

r r J 


4 


(¥) 


(3 - 4i/) 
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SGpj _ _1_ 

dik ” 4?r \k(X + 2/i) 

dFjj _ 1 1 ;| 

47rr 2 (1 — t/) 


[(*?)*> -<*){t + t}] 

4 yiyjykyini m/j^k jjkymni 


+ 


r* 

"ZyjWSH 


t f-\ ( MijykVini c _ , 2yiy k nj t _ \ 7 

/h 7 ?) ( r 2 bij n k + ^2 ~ ) 72(7?) 


“ /sfa) 


- £ (rfe) [(^) «*■ - '-"'l - (SL + SO. + f«2i) ft) 


d^k 
fi(v) = 2 
/2 M = 1 - 2j/ 


= 1 — 2*/ 
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APPENDIX B.2 - 


Kernels for Steady Incompressible Thermoviscous Flow 


4tt/i 


ViVj 


— Sijln r 


F l3 = 


1 

2nr 


^yiVjVknk 

r 3 


j 1 k Vi ^ ^ikVj Vk 

dxk 4njir r r r 


^yiVjVk 

r 3 


GM = 2^ [/n r ] 



VkTlk 

r 


dGee _ 1 yk 

dxk 27r kr r 


yi = Xi - fc 


r 2 = y*yt 
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APPENDIX B.3 - 

Kernels for Unsteady Incompressible Viscous Flow 


Gait-X't) 


4tt M * ij '{ 


2 




Fiji! -X,t)= - e -«*/<} + ^~-{si(t]) - H(t )} + - e^*} 


-*,<) = 1 


3x fc 


47 r/ir 


^*iM) + ^<..(1)} - - »,M) 




where 


Vi=ti- Xi r 2 = yij/j 

*> = (Syr7= c = p/p 

si(v)=%{ l-e"^) 

£i(*) = JT^ u . 


Then, 

GgK-^sG^-X.nAr) for n = l 

G^-X) = G^-X >n Ar)-G fi (e-X,(n-l)Ar) for n>l 


with similar relationships for - *) and Sf K - *)• 
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APPENDIX B.4 - 

Kernels for Steady Convective Incompressible Viscous Flow 


Gpj = ~2 w (r) 

Fij = fi Uk + G PJ' n< + P U kGijTl k 


W) e -fi K < a) _ 1 ti) tt _ £ ( Ei\ tt , L Mi) ft 

u 2 J u\u) dxj u \ u ) dxi u \ u ) dx t 



/ c \ (Uj\ d 2 <p /c\ fSijUA d 2 4> ' 

\uJ V £/ J dxidx k + Vf/J \ u ) dx e dx k . 


where 


yi = Xi- ii, 


r 2 = ym 


P = U k y k /2c 


U 2 = UiUi 


a — U r/2c 

<j> = - ln(a ) - e~ p K 0 (a) 


dp_ 

dxi 



+ 



e- p K 1 (a) + 



e p K 0 (a) 


ay 

dxidxj 




8 2 <(> 
dxjdxk 


e ^K 0 (a) 
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APPENDIX B.5 - 

Kernels for Steady Convective Compressible Thermoviscous Flow 


47r(A + 2p) 


1 6 ij e u "»' i '>Ko 6ijU k y k - U iVj - U jVi ) 


+ 4rfi { 6 ij eUkVk/2 ‘' K ° (1^) + ^ii U kVk ~ Ujyj U jVi ) 


2v 

Ur 


— e 




Gip — G\ p + G 2 p 


Hjc-U) 
lp _ 2ir Po U 2 [ 

H(U - c) c 


v / (WW - ^ 


+ 


2 p 0 I/ 2 v 

H{c-U) 

2p ~~ 2* Po U 2 [ 


-U x H{U k y k - vr) 


GL = - 


iu. In v'(W + »V - * U^, B } )] 


+ ^jTT- ^HiUkyk-vr) 

2p 0 U 2 v 



1 

2irp 0 U 2 


H(c - U) 
:pp ~ 2tt 


Ui^y^Kp + ±-{U iyk y k - U 2 yi ) r e u ^r Ko{m 

\2r) J Ur JUr/2r) 

; >” V , («kS») s + « 2 f ! - V»»)* +•’>■’] -\ H( - U ~ c )z H ( Ukm ~ '" r) 


f? = (A + 2/i)/po, ^ = p/Po , K = 


t> 2 = |c 2 — f/ 2 | 
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